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SUMMARY 


The Rotorcraft Dynamics Division, Aeroflightdynamics Directorate, U.S. Army Avia- 
tion Research and Technology Activity (AVSCOM) has developed the General Rotorcraft 
Aeromechanical Stability Program (GRASP) to calculate aeroelastic stability for rotorcraft 
in hovering flight, vertical flight, and ground contact conditions. In this report, GRASP is 
described in terms of its capabilities and the philosophy behind its modeling. The equa- 
tions of motion that govern the physical system are described, as well as the analytical 
approximations used to derive the equations. These equations include the kinematical 
equation, the element equations, and the constraint equations. In addition, the solution 
procedures used by GRASP are described. 

GRASP is capable of treating the nonlinear static and linearized dynamic behavior 
of structures represented by arbitrary collections of rigid-body and beam elements. These 
elements may be connected in an arbitrary fashion, and are permitted to have large rel- 
ative motions. The main limitation of this analysis is that periodic coefficient effects are 
not treated, restricting rotorcraft flight conditions to hover, axial flight, and ground con- 
tact. Instead of following the methods employed in other rotorcraft programs, GRASP is 
designed to be a hybrid of the finite-element method and the multibody methods used in 
spacecraft analyses. GRASP differs from traditional finite-element, programs by allowing 
multiple levels of substructures in which the substructures can move and/or rotate relative 
to others with no small-angle approximations. This capability facilitates the modeling of 
rotorcraft structures, including the rotating/nonrotating interface and the details of the 
bla.de/root kinematics for various rotor types. GRASP differs from traditional multibody 
programs by considering aeroelastic effects, including inflow dynamics (simple unsteady 
aerodynamics) and nonlinear aerodynamic coefficients. 
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1. INTRODUCTION 


Previous helicopter aeroelastic stability programs have suffered from significant re- 
strictions. The General Rotorcraft Aeromechanical Stability Program has been developed 
using a modern approach which overcomes these limitations. 

1.1. Background 

In early efforts made to calculate the aeroelastic stability of hingeless helicopter ro- 
tor blades, it was common practice to make use of simple physical models ( e.g ., spring- 
restrained, centrally-hinged, rigid blades (ref. 1)). Later work treated configurations that 
were somewhat more complex, and included models of elastic blades (ref. 2), body degrees 
of freedom, and inflow dynamics (ref. 3). These simple approaches to rotorcraft aeroelastic 
stability calculations have been very valuable for gaining physical insight into many com- 
plicated phenomena (e.g., coupled rotor-fuselage stability). They all are, however, based 
on a single physical model, and therefore are of limited value when more realistic rotorcraft 
configurations must be analyzed. 

Because of the complex couplings inherent in a bent and twisted beam, the calculation 
of aeroelastic stability is particularly important in the analysis of rotor blades having 
cantilever root boundary conditions (e.g., hingeless and bearingless rotors). In bearingless 
rotors, the blade/root kinematics demand a great deal of modeling flexibility because 
individual blade designs tend to have widely varying configurations. The FLAIR program 
(refs. 4, 5, and 6) is able to perform this type of aeroelastic stability calculation, but is 
limited to a configuration that has a rigid blade, a uniform flexbeam, linear aerodynamics, 
static induced velocity, and several different blade/root configurations. While FLAIR is 
currently being used in the rotorcraft community, it lacks the flexibility and generality 
necessary for it to be considered general-purpose analysis. 

For analysis of problems involving complete rotorcraft, there exist large helicopter sim- 
ulation programs such as C-81 (ref. 7) and G400 (ref. 8). These programs were designed 
primarily for time-history analysis of rotorcraft behavior in forward flight rather than for 
aeromechanical stability. Despite their generality and complexity, these programs have 
limitations (primarily related to aerodynamics) which are pointed out by Johnson (ref. 9) 
in his discussion of these and other large rotorcraft programs. While the CAMRAD pro- 
gram overcomes m.my of these limitations, all of these programs (including CAMRAD) are 
restricted to a fixed number of physical models, and lack the modeling flexibility needed to 
deal with a wide variety of blade/root geometries. Many of these programs rely on results, 
such as a set of inodes, from other programs. This approach may present an assortment of 
modeling difficulties, especially for bearingless rotor blades. In particular, the mathemati- 
cal and physical consistency of a combined approach is seldom examined, and the physical 
bases of the individual programs are likely to not be consistent. Furthermore, in stability 
analyses a nonlinear static equilibrium solution is needed about which to linearize — an 
important consideration which most of the earlier simulation programs do not address. 
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Therefore, it is important that a code be developed in which blade structural dynam- 
ics, isolated blade stability, and isolated rotor stability, as well as coupled rotor/airframe 
stability, can all be treated under a consistent set of physical assumptions. 

Dynamic coupling programs, such as DYSCO (ref. 10), which have a high degree 
of generality, allow coupling of discrete component models and/or modal representations 
of flexible structures. While DYSCO has a very powerful, executive- driven system, it 
currently cannot treat the aeroelastic behavior of bearingless rotor systems undergoing 
geometrically nonlinear deformation. The problem is that it lacks a sufficiently general 
element in its element library. 

Several recent implementations that apply the finite-element method to rotorcraft 
problems (refs. 11, 12, and 13) are not able to overcome these limitations because their 
physical models are limited to a single configuration. Simply breaking a rotating beam 
into a number of finite elements yields nothing more than a discretized rotating beam. 
This approach does not meet the requirement that the beam be coupled with an airframe, 
or model blade/root kinematics of an arbitrary configuration. The classical finite-element 
method is based on the breaking up of a single structure (i.e., a beam, plate, or shell) 
into an arbitrary number of elements and expanding the appropriate field variables into 
polynomial shape functions. This approach by itself also lacks the flexibility to deal with 
truly arbitrary rotorcraft configurations because a helicopter is a system of structural 
components, some of which may be rotating and/or translating relative to one another. 
Because of this, rotorcraft are actually more akin to the multibody systems (refs. 14 and 
15) encountered in spacecraft problems. Unfortunately, few multibody programs possess 
the capability to deal with flexible components, and none have the capability to deal with 
aeroelastic phenomena since they were developed primarily for spacecraft applications. 

All previous attempts at modeling rotorcraft problems have incorporated certain re- 
strictions that are undesirable in a truly general-purpose program. General-purpose codes 
that are currently under development, or will be developed in the future, should over- 
come the major shortcomings of existing aeroelastic analyses. Consider, for example, the 
following typical restrictions: 

The first is a restriction to linear, small-displacement approximations of beam elastic 
deformation. This restriction is unacceptable in a general-purpose rotorcraft program 
because the rotor blade aeroelastic problem, especially for hingeless and bearingless rotor 
blades, has been conclusively shown to be a nonlinear problem. A consistent approach 
based on nonlinear kinematics is required for these configurations. 

The second is a restriction to elastic blade models with ordering schemes, second- 
degree nonlinearity, or “moderate” rotations. These approximations are undesirable be- 
cause the governing equations often have to be augmented with certain higher-order terms 
if the values of certain structural properties are not within some nominal range (Rosen 
and Friedmann (ref. 16)). Therefore, in a general-purpose analysis, the higher-order terms 
must be present. Ordering schemes, while still a valuable tool when used in special-purpose 
codes and codes where accuracy is a secondary consideration, are neither necessary nor 
desirable in a general-purpose context. Furthermore, a bearingless-rotor flexbeam must 
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undergo deformation-induced rotations of the order of the collective pitch angle — a ro- 
tation too large to be classified as “moderate.” Thus, bearingless rotor problems demand 
a large deflection analysis without artificial restrictions on rotations due to deformation, 
the degree of nonlinearity, or the values of blade properties. 

The third restriction is to a fixed number (usually one) of configurations ( e.g isolated 
hingeless blade or coupled bearingless rotor and body or a single blade/root configuration). 
This restriction is unacceptable in a general-purpose code because the intent of such a code 
is to analyze different types of configurations with a single, consistent set of assumptions. 
Such a code should be able to treat all currently known blade/root mechanisms and, at the 
same time, model configurations that do not yet exist. It should be possible to construct a 
new configuration with simple building blocks and with no artificial limitations on the pro- 
cess. For maximum flexibility in treating these different configurations, the finite-element 
method is the preferred approach. Moreover, the existence of many different blade/hub 
configurations for helicopters requires a capability to analyze arbitrary configurations of 
structures, parts of which may be rotating. Thus, the code should employ the multibody 
philosophy. 


1.2. Approach 

To overcome the aforementioned limitations of the existing methods of aeroelastic 
stability analysis, the General Rotorcraft Aeromechanical Stability Program has been de- 
veloped. GRASP combines the finite- element and multibody approaches, and incorporates 
multiple levels of substructures to provide a powerful tool for rotorcraft analysis. The de- 
sign of GRASP is based on the concept of a collection of flexible and rigid bodies connected 
in an arbitrary manner. Libraries of elements, constraints, and solution algorithms appro- 
priate for the helicopter aeroelastic stability problem were designed and built into the 
program. 

The element library promotes the modeling of the blades as beams; construction of 
arbitrary mechanisms to treat blade/root. kinematics with beam elements and rigid bodies; 
treatment of the fuselage as either a rigid body, a collection of beam elements, or a modal 
representation obtained from some other source; and treatment of both static and dynamic 
induced inflow by means of blade-element /momentum theory. The constraint library allows 
arbitrary connections between elements, includes constraints that allow for compliance in 
the constrained relative motion between elements, and includes constraints that allow 
the connection of rotating and nonrotating substructures. None of the constraints in the 
library use any kinematical approximations, such as small-angle assumptions. The solution 
procedures include nonlinear static equilibrium and linearized stability about equilibrium, 
both presently limited to the hovering flight condition. 

It should be noted that these physical modeling assumptions and solution procedures, 
while adequate for aeromechanical stability analysis in axial flight and ground contact, 
are not adequate for a comprehensive rotorcraft dynamic analysis as defined by Johnson 
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(ref. 9). The analysis methodology used in GRASP, although a viable approach for ap- 
plication to nonlinear dynamics in forward flight, would require considerable effort to be 
implemented in GRASP. 

Several very desirable, but not required, features of a general-purpose code, have 
been incorporated in GRASP. 1) The accuracy of the analysis may be increased without 
having to add more elements. The aeroelastic beam finite element developed specifically 
for GRASP uses a variable-order (or p-version) approach, which is based on high-order, 
orthonormal, polynomial displacement functions (refs. 17 and 18). 2) As much as possible, 
the equations of motion are formed by the program internally, minimizing the possibility 
of human error in the equations. 3) The user interface is capable of handling a general 
problem without having to be supplied with the form of the equations of motion or even 
the number of degrees of freedom. 4) Both large and small problems can be modeled 
with the same code. Thus, the number of degrees of freedom is not fixed a prion. This 
feature not only requires a great deal of flexibility in assembling the system equations of 
motion, but also requires that data be structured and managed in core with a flexibility 
not, inherent in FORTRAN (ref. 19). 
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2. SOLUTION APPROACH 


GRASP is specifically designed to provide a tool for determining the equilibrium 
deflections and aeroelastic stability of arbitrary rotorcraft configurations in hover or vertical 
flight. A GRASP rotorcraft model is considered to be an aeroelastic system consisting of 
a structural system, portions of which may be rotating relative to one another, and a 
moving air mass with which the structure interacts. All parts of the model may be subject 
to forces and externally applied constraints. The position of any point on the structure or 
the air velocity at any point in the flow field relative to an inertial frame of reference may 
be determined by solving a system of partial differential and boundary value equations. 
These equations are obtained from the laws of fluid and structural mechanics, and from 
the constitutive properties of the materials in the structure and the air. 

In vertical flight, hover, or ground contact a rotorcraft can assume a steady-state 
equilibrium configuration when the airflow, gravity, and the rotor spin axis are aligned; 
and when the angular velocity of the rotationally isotropic rotor is constant. In this 
restricted case where the structure is not subject to time-varying forces, it is possible to 
eliminate explicit time dependence from the equations. This steady-state equilibrium can 
be considered to be static when contrasted to the more general periodic equilibrium found 
in forward flight problems. The steady-state equilibrium configuration is characterized by 
a time-invariant deformation in the nonrotating portions of the rotorcraft, a steady flow of 
air through the rotor disk, and time-invariant deformations of the rotor blades with respect 
to a rotating reference frame. The steady-state solution then calculates the equilibrium 
values of all of the model generalized coordinates and generalized forces. 

The equations of motion for the continuous-structure portions of the structure are 
discretized by means of variable-order, finite-element shape functions. The equations for 
the structure then become a system of nonlinear, ordinary differential equations. It is 
possible, as indicated above, to eliminate all explicit dependence on time from the equations 
for the restricted case of axial flight or ground contact. A linearized system of equations 
may then be calculated by taking small perturbations about the static equilibrium state. 
The stability problem is defined, therefore, by a second-order system of linear equations 
with constant coefficients. 

For infinitesimally small perturbations about a previously-calculated, steady-state 
configuration, the dynamic motion of the rotorcraft can be represented as a linear combi- 
nation of complex eigensolutions. Since the aeroelastic stability of the rotorcraft can be 
determined directly from the eigenvalues, the primary objective of GRASP can be satisfied 
by computing these eigensolutions. The frequency and damping information in the eigen- 
values and the modal information in the eigenvectors, which can also be obtained from the 
eigensolutions, facilitate the user’s understanding of the dynamics of the rotorcraft. 

The eigensolution provides the complex eigenvalues and eigenvectors for all model 
degrees of freedom associated with the equations of motion Mq + Cq + Kq ~ 0 which have 
been linearized about a steady-state deformation. These equations are often referred to 
as being “asymmetric” because of the nonsymmetry due to aerodynamics contributions 
to the coefficient matrices C and K . The coefficient matrix A/, which is both symmetric 
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and positive-definite, contains contributions from the mass of the structural model and 
from the “apparent mass” of the air. The coefficient matrix C contains contributions from 
structural and aerodynamic damping and inertial forces. The coefficient matrix K contains 
contributions from structural stiffness and effective stiffness from aerodynamic and inertial 
forces. Like the steady-state solution, this solution requires that the model correspond to 
a physical system which is not subject to time- varying forces. Currently, the asymmetric 
eigensolution must be computed by using the steady-state solution obtained for an identical 
model. This solution procedure prohibits one, for example, from obtaining the steady-state 
deformations of an isolated blade, then applying that solution to a coupled, rotor/fuselage 
configuration. 
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3. MODELING APPROACH 


In order to form a mathematical representation of a structure that may contain bodies 
which are experiencing large kinematic motions relative to one another, it is necessary to be 
able to write the full, nonlinear equations of motion for the structure. The fundamentals of 
the approach used in GRASP to derive these equations are adapted from methods that were 
originally developed for spacecraft applications (ref. 20). For the types of structures that 
GRASP is designed to represent, additional emphasis has been placed on using multiple 
levels of substructures to model a structure. 

The first step in modeling a structure in this manner is to decompose the structure 
(called the parent) into a set of subordinate substructures (called children), each of which in 
turn may also be decomposed into a set of child substructures. This decomposition process 
continues until every substructure has been decomposed into simple structural elements. 
The lowest-level substructures ( i.e . those with no children) are called elements. The result 
of this method of modeling a structure is a hierarchically-ordered set (tree) of substructures 
(fig. 1) that has the complete structure at the root and elements at each of the leaves. 
Under this modeling scheme, a parent substructure may have any number (including zero) 
of child substructures but only one parent substructure. The only substructure without a 
parent is the complete structure, which is at the root of the tree. 



Model-type subsystem: 1 
System -type subsystems: 2, 3 
Element type subsystems: 4, 5, 6, 7 


Figure 1. Hierarchical substructure tree. 

The hierarchical model representation implemented in GRASP allows great generality 
in the types of configurations that can be analyzed, and permits essentially arbitrary kine- 
matic motions of components relative to one another. This general framework, along with 
a software design that emphasizes the use of libraries for constraints, elements, solutions, 
and so on, means that the capabilities and limitations of the program are those associated 
with the members of the libraries, not with the program in general. 
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3.1. Subsystems 

In GRASP, substructures are abstracted into subsystems. Each substructure is then 
represented by a subsystem, which may be classified according to its position in the hi- 
erarchy (fig. 2). The subsystem representing the complete structure (or model) is called 
a model-type subsystem. Substructures having no children (elements) are represented by 
element-type subsystems. The remaining subsystems, those having a parent and at least 
one child, are represented by system-type subsystems. To represent the substructures that 
make up the model, subsystems serve several functions. First, they contain the complete 
definitions of the substructures that they represent. Second, they are repositories for the 
generalized coordinates, generalized forces, and dynamic matrices associated with the sub- 
structures. Finally, they serve as the basic units of the hierarchical organization, which 
is an integral part of the computational process of transforming the parent generalized 
coordinates to the child generalized coordinates, and transforming the child generalized 
forces to the parent generalized forces. 



Figure 2. Hierarchical subsystem tree. 


Subsystems, in general, may contain the following: a frame of reference, a set of nodes, 
a set of generalized coordinates, and a set of constraints. Each of these entities performs 
a different function within the subsystem, and will be described in the following sections. 
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3.1.1. Frames of Reference 


Every subsystem in a GRASP model (with the sole exception of the air mass element) 
has a frame of reference associated with it. The frame of reference is not associated with 
any material point on the substructure, but instead serves as the “point of view” for the 
subsystem. As such, it establishes the coordinate system for that subsystem. The initial 
position and orientation of a reference frame may be selected to define a coordinate system 
that is natural for the subsystem (e.g., a hub-centered frame of reference might be selected 
for a subsystem that contains a helicopter rotor). 

Since reference frames are not physically connected to any structure, but rather are 
allowed to move freely, six degrees of freedom are associated with each frame. These 
degrees of freedom define the position and orientation of the frame of reference for the 
current subsystem relative to the reference frame for its parent subsystem. 

In addition to serving as a reference for the subsystem, the frame of reference may be 
used to model the discrete motions of the substructure. This can often lead to significant 
simplifications in the equations of motion for subordinate subsystems. For example, if a 
reference frame is attached to the root of a rotating beam and used to model the rotational 
motion of the beam, the equations of motion of the beam itself need not explicitly include 
the rotational motion. 

Since Newton’s laws hold only in an inertial reference frame, the model-type subsystem 
at the root of the tree is defined to be fixed in an inertial frame of reference. Therefore, 
while a model-type subsystem does have a frame of reference, that reference frame has no 
degrees of freedom associated with it since it must be inertial. As a result, the motions of 
every part of the system can be related to an inertial frame of reference. 

3.1.2. Nodes 

Nodes are used by GRASP to introduce degrees of freedom into a model. In general, 
the degrees of freedom introduced by a node may be any generalized coordinates that can 
be associated with a physically identifiable property of the structure. For example, the set 
of degrees of freedom for a node could be defined to be the three rigid-body translations 
and the three rigid-body rotations of a point on a structure. Alternatively, there could be 
a node whose degrees of freedom are defined to be modal coordinates. 

Currently, two different types of nodes are used by GRASP: structural nodes and air 
nodes. The structural nodes provide the measures for the local displacement and rotation 
of a structure. They move with the deformation of the structure and may be conceptualized 
as massless, infinitesimal, rigid bodies that are physically attached to the structure. The 
air nodes define the induced inflow velocity field through a helicopter rotor. The degrees of 
freedom for the air node are measures of the velocity distributions around the rotor disk. 
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3.1.3. Constraints 

The constraints act as a sort of “glue” that holds a model together. Constraints 
are used to model both physical constraints ( e.g ., pins, gimbals, and clamps), and to 
eliminate the dependent degrees of freedom that have been introduced into the model. An 
example of a physical constraint would be the clamped boundary condition at one end of 
a cantilever beam. That end of the beam is modeled by constraining the node at one end 
of the beam to have no translational or rotational motion. Now consider two frames of 
reference that are defined to move as if they are rigidly connected to one another. For this 
system, there are twelve degrees of freedom (six for each frame), but only six of them are 
independent. Therefore, a constraint must be defined to remove the dependent degrees of 
freedom. In general, the set of constraints for a subsystem must be sufficient to reduce the 
total number of degrees of freedom to only the independent degrees of freedom for that 
subsystem. Similarly, for the complete model, all dependent degrees of freedom must be 
eliminated. 

All of the constraints implemented in GRASP are based on purely kinematical rela- 
tionships. There are no restrictions to small or moderate displacements or rotations in any 
of the constraint equations. However, it is necessary to avoid the singularity that occurs 
for deformation-induced rotations of 180° . This singularity arises as a result of using 
finite- rotational kinematics that are based on Rodrigues parameters (ref. 21). 

The constraints in GRASP are implemented at two levels: the program level and 
the user level. The constraint “primitives” are found at the program level. These simple 
constraints provide a basic set of connections among generalized coordinates, frames, and 
nodes. At the user level, these primitive constraints are combined to provide the user with 
physically meaningful constraints between structural elements. For example, the rigid- 
body mass connectivity constraint, which is used to attach a rigid-body mass element to 
a structure, is a combination of a primitive constraint between frames and a primitive 
constraint between nodes. 

In order to provide a full set of constraints, the constraint library in GRASP includes 
several different classes of constraints. These include constraints between two frames, con- 
straints between tw r o nodes, constraints between generalized coordinates, and constraints 
between a, frame and a node. 


3.2. Elements 

Elements are subsystems that have no children. In addition to frame and nodal 
degrees of freedom, they may also have additional, non-nodal generalized coordinates. 
Computationally, the elements are the primary source of virtual work in the structure. 
For steady-state problems, the elements return the generalized forces associated with a 
given set of generalized displacements. For perturbation problems, the elements return 
the element coefficient matrices. These matrices are determined from the perturbations 
in generalized forces resulting from perturbations in the generalized coordinates and their 
time derivatives. 
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3.2.1. Rigid-body Mass 

The rigid-body mass element represents a rigid body that is subject only to inertial 
and gravitational forces. It has a single structural node that is located at the mass center, 
and its axes are aligned with the principal axes of the body. The frame of reference for the 
rigid-body mass element coincides with the nodal coordinates in their undeformed state. 

3.2.2. Air Mass 

The air mass element models the momentum air flow through an axisymmetric rotor 
disk. The degrees of freedom associated with this element are introduced through a single 
air node. Since the air mass element is defined to be fixed in inertial space, the frame 
degrees of freedom are suppressed. For steady-state problems, the residuals corresponding 
to the uniform inflow velocity and the radial velocity gradient are calculated from momen- 
tum considerations (ref. 22). For the asymmetric eigenproblem, only the momentum terms 
(ref. 23) involving uniform and first- harmonic, cyclic perturbations of the inflow velocity 
contribute to the element coefficient matrices. 

3.2.3. Aeroelastic Beam 

The aeroelastic beam element represents a slender, nonuniform beam (without shear 
deformation) that is subject to elastic, inertial, gravitational, and aerodynamic forces. 
The primary assumption in the derivation of the element equations (ref. 24) is that strains 
remain small relative to unity. There are no small-angle approximations made and all 
kinematically nonlinear effects are included. One current limitation is that orientation 
angles (ref. 21) (of type body three: 1 2 3) are used in the description of finite rotation 
inside the beam element. Thus, rotations due to the deformation of beam elements may 
not exceed 90° . 

The aeroelastic beam element degrees of freedom come from a frame of reference that 
coincides with the root of the element in its undeformed state, structural nodes at the 
root and tip, an air node, and a set of internal degrees of freedom. The internal degrees 
of freedom result from the higher-order polynomials that may be used to increase the 
accuracy of the beam deformation calculations. When no internal degrees of freedom are 
specified, the aeroelastic beam is an Euler-Bernoulli beam in which the axial and torsional 
deformations in excess of a built-in pretwist are represented by linear polynomials, while 
the bending deflect : ons are represented by cubic polynomials. The method of adding 
internal degrees of freedom to improve the accuracy of an element is more convenient than 
adding elements, and is also more efficient (ref. 17) given the same number of degrees of 
freedom. Internal degrees of freedom may be added selectively to reflect the dynamics of 
the clement. For example, if a beam is very stiff in bending and extension but soft in 
torsion, additional torsional degrees of freedom may be added without having to include 
any more bending or extensional degrees of freedom. 

The aerodynamic forces on the beam element are calculated from quasi-steady strip 
theory using lift, drag, and moment coefficients that are piecewise continuous functions of 
the angle of attack. Spanwise scale factors for the lift, drag, and moment may be specified 
to allow for tip loss and other similar effects. The chord width, the pitch angle of the 
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zero-lift-line, and the offset of the aerodynamic center from the elastic axis may also vary 
over the length of the element. The aeroelastic beam element also calculates the blade- 
element contributions to the induced velocity, which are combined with the momentum 
contributions from the air mass element elsewhere. 
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4. SOLUTION METHODS 


The solutions currently implemented in GRASP allow the user to calculate the steady- 
state deformations of a structure under load, and then to solve for the eigenvalues and 
eigenvectors of the deformed structure. In order to obtain a valid eigensolution, the steady- 
state deformations that are used must be such that the structure is in equilibrium. 

4.1. Steady-State Solution 

The equations for the steady-state equilibrium of the model are a set of nonlinear, 
algebraic equations of the form 

Qi = /(?!, i = (4.1 l) 

where the Qi are the generalized forces (residuals), the qi are the generalized coordinates, 
and N is the number of system degrees of freedom. These equations are generated in- 
ternally by GRASP at the element level, and automatically assembled by the constraints, 
which combine the contributions from the finite elements into the final set of equations. The 
solution to this set of equations is obtained through the use of the Levenberg-Marquardt 
algorithm. This algorithm minimizes the sum of the squares of the residuals from the 
steady-state equations. The implementation in GRASP uses the IMSL (ref. 25) subrou- 
tine ZXSSQ. 

For problems involving the aeroelastic beam element with internal degrees of freedom, 
the solution algorithm is used at two levels. First, it is used in an outer iteration loop to 
arrive at a solution to the steady-state equations for the complete model (which excludes 
the aeroelastic beam internal degrees of freedom). In addition, it is used in a separate, 
inner iteration loop to calculate the internal degrees of freedom for each aeroelastic beam 
element. A full inner solution for each aeroelastic beam is calculated for each iteration of 
the outer solution. 

In order to arrive at a steady-state solution, the residual forces on the system must 
be calculated, given a deformation state. The algorithm that is used to calculate the 
residuals for the top-level subsystem in the hierarchical organization of the model is based 
on a full-order tree traversal (fig. 3). When traversing down the tree (away from the root 
subsystem), the st:~i,e vector for each child subsystem is calculated from that of its parent. 
Also, the inertial motion of the child subsystem reference frame is calculated from that of 
the parent. Upon reaching an element, the state vector for that element and the inertial 
motion for the element frame are used to calculate the element residuals. Traversing 
back up the tree (towards the root subsystem) the residuals from each child subsystem 
are transformed into its parent subsystem and added to the parent residuals. When the 
traversal is complete, the residuals corresponding to each generalized coordinate in the root 
subsystem are known. The complementary processes of calculating the state vectors and 
assembling the residual vectors are accomplished by using the constraints, which define 
the relationships among the degrees of freedom in the parent and child subsystems. 
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Calculate generalized coordinates (1, 2, 3, 5, 7, 10) 

Assemble generalized forces (4, 6, 8, 9, 11, 12) 

Forces calculated in subsystems 4, 5, 6, 7 

Figure 3. Steady-state solution full-order traversal. 

The solution methods available in the current version of GRASP are restricted in that 
the same model must be used for both the steady-state and asymmetric eigenproblem solu- 
tions. This creates a problem for the steady-state solution algorithm when a configuration 
contains unconstrained degrees of freedom. This can occur when a. model having both ro- 
tating and nonrotating components is being analyzed. For such a configuration, the cyclic 
degrees of freedom generated by the rotating constraints are unconstrained. It can also 
occur in airborne configurations, which suffer from the same problem because their body 
degrees of freedom are unconstrained. To alleviate this problem, GRASP currently marks 
these unconstrained degrees of freedom during the building of the model, and eliminates 
them from the state vector used in the minimization algorithm. 


4.2. Asymmetric Eigenproblem Solution 

The system equations for the asymmetric eigenproblem can be expressed in the famil- 
iar form 

M + C'q + Kq = 0 (4.2-1) 

where the q ' s are infinitesimal perturbations of the generalized coordinates. The algorithm 
used to assemble the coefficient matrices for the root subsystem is very similar to that 
used to calculate the steady-state residuals in that it also is based on a full-order tree 
traversal (fig. 4). However, while traversing down the tree, no state vector calculations 
are required. Upon reaching an element, the coefficient matrices for that element are 
calculated. During the traversal back up the tree, the constraints are used to assemble the 
child subsystem matrices into the parent matrices. At the conclusion of the traversal, the 
coefficient matrices for the model subsystem are complete. 
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Assemble subsystem matrices (4, 6, 8, 9, 11, 12} 
Calculate element matrices in subsystems 4, 5, 6, 7 


Figure 4. Eigensolution full-order traversal. 


The solution of this set of equations is begun by factoring matrix M using the Cholesky 
decomosition algorithm. The GRASP implementation uses subroutine LUCECP from the 
IMSL (ref. 25) library. M then becomes 


M - LL J 


(4.2-2) 


Introducing the transformation 


L T q 


(4.2 3) 


the mass matrix M can be reduced to an identity matrix and the system equations can be 
written as 

Ai + L~ 1 CL~ t z 4- L~ l K L~ T z = 0 (4.2-4) 

Writing this system of equations in first-order form 


A 

0 ' 

. 

0 

A 1 

0 

A 

y = 


1 

0 
7 

1 


where 



Time may be eliminated by the introduction of 



(4.2-5) 


(4.2-6) 


(4.2-7) 
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which allows the extraction of eigenvalues and eigenvectors directly from the matrix on the 
right-hand side of equation (4.2-5). The dynamic matrix is balanced, converted to Hes- 
senberg form, and then the QR algorithm is used to obtain the eigensolution. Finally, the 
eigenvectors are transformed back to the original coordinate system via the transformation 

q* = L T z * ( 4 - 2 - 8 ) 

GRASP uses subroutine RG from the NASA/Ames Cray library to calculate the 
eigensolutions . 
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5. COORDINATE SYSTEMS 


In GRASP, many different coordinate systems are used to mathematically describe 
the physical structure. To differentiate among them, each coordinate system is identified in 
its undeformed state by a capital letter ( e.g ., A), Depending on the context, an identifier 
may refer either to the coordinate system itself or to a point located at the origin of the 
coordinate system. The addition of a prime or a double-prime to the identifier indicates 
that the designated coordinate system either is in a state of static equilibrium (e.p., A 1 ) or 
is in a dynamically perturbed state A”). With these multiple coordinate systems, it 
is often desirable to use several types of mathematical notation when deriving and writing 
equations. Not only can the form of the equations be simplified, but also they can be 
made more readable. This section is intended as an introduction to the notation used in 
the sections where the equations are actually derived. 

5.1. Vectors 

Vectors play an important role in coordinate system mathematics. Associated with the 
orthogonal axes emanating from the origin of every coordinate system is a set of dextral 
unit vectors. These unit vectors are called the base or basis vectors of the coordinate 
system. In addition, vectors are used to define variables such as position, velocity, and 
acceleration. Three types of notation are used in writing vector expressions and operations: 
vector-dyadic notation, index notation, and vector notation. 

5.1,1. Vector- Dyadic Notation 

All vectors and dyadics used in GRASP are underlined V^), and all unit vectors 

are identified by a circumflex. The difference between a vector and a dyadic should always 

be clear from the context of its usage. For a coordinate system A, the basis vectors are 
-A 

written as b ± , where i == 1,2,3. Any unit vector other than a basis vector is denoted by e, 
and may appear either with or without superscripts. 

When kineinatical quantities have coordinate systems associated with them, the rela- 
tionship is defined by using the appropriate superscripts. For example, 

R BA = position of the origin of coordinate system B 

with respect to the origin of coordinate system A 

V s A = velocity of the origin of coordinate system B 
with respect to coordinate system A 

A ba = acceleration of the origin of coordinate system B 
with respect to coordinate system A 

Q BA = angular velocity of coordinate system B 
with respect to coordinate system A 
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Forces and moments are significant in their point of application as well as their source. 
The notation adopted herein is 


For example, 
relationship 


F a =force at A 
M_ a ^moment at A 

a force and moment at A contribute to a moment at D according to the 

M b =M a + R ab X A a (5.1 .1—1 ) 


5.1.2. Index Notation for Vectors 
A vector V in the A basis may always be expressed as 


V = V Ai b? (5. 1.2-1) 

where the summation convention adopted is that repeated indices are always summed 
over their range. Unless otherwise specified, Latin indices assume the values 1,2,3; Greek 
indices assume the values 1,2. The subscript A in V Ai indicates that the measure numbers 
VjLi are defined in the A basis. 


Two symbols frequently encountered in vector operations that 
the Kronecker delta b x j and the Levi-Civita epsilon tijk where 

use index notation are 


6» = ' 

[ o i^j 
\l i=j 

(5. 1.2-2) 


' 0 

any index repeated 


^ ijk ^ 

+1 

cyclic permutation 

(5. 1.2-3) 


-1 

V. 

acyclic permutation 



The Kronecker delta consists of the components of the identity tensor in a Cartesian 
coordinate system, while the Levi-Civita epsilon consists of components of the permutation 
tensor in a Cartesian coordinate system. Some useful identities regarding both of these 
symbols may be found in reference 26. 


5.1.3. Matrix Notation for Vectors 


Using index notation, a vector V may be expressed in the A basis as shown in equa- 
tion (5. 1.2-1). Since the basis is identified by the subscript A, the measure numbers 
themselves may be viewed as a complete description of the vector. Thus, the column 
matrix Va can be defined to be 

f V A1 ) 


V A 



(5. 1.3-1) 
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as an alternate way of expressing the vector ]£. The dot product 17 • F may then be written 
as 


'A 1 

cJva - l^i u A2 u A3 1 <( r A2 


V'. 


A3 


(5. 1.3-2) 


The cross product of two vectors U_ and V_ may be written as 


U x V = U Ai b* x V Ai b? 


— e ijk U Aj V Ak b, 


(5. 1.3-3) 


u Aij v Aj b? 


This equation implies that the measure numbers of the cross product in the A basis are 
simply the elements of the matrix product U A V A where 


For example, 


( )ij ~ C ijk{ )fc 


U A = 


0 

V 43 

-U A 2 


~U A 3 
0 

Uai 


U A 2 

-Uai 

0 


(5. 1.3-4) 


(5. 1.3-5) 


There are also several useful identities that can be derived for two column matrices a 
and b 

-T 

a = — a 


ah = — ha 


a T b = — b T a = (ab) T 
ab = — a 1 bA 4- ba 7 


(5. 1.3-6) 


ab —ba T — ab T 


ab 


ab =ba 
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5.2. Finite Rotations 


In many kinematic analyses, rotations are assumed to be either infinitesimal or moder- 
ate in size. These assumptions allow certain simplifications in the kinematical relationships, 
but constrain the range of applicability of the analysis. In GRASP, no such assumptions 
are made and all rotations are assumed to be of arbitrary size (finite). Finite rotations are 
expressed in four ways in GRASP: 

(1) direction cosines, 

(2) Euler rotations, 

(3) Tait-Bryan orientation angles, and 

(4) Euler-Rodrigues parameters. Internally, GRASP expresses all finite rotations in terms 
of direction cosine matrices. For the convenience of the user, any of the other three 
methods may be used to specify the input to GRASP. Since there are significant 
differences in the algorithms used to compute the direction cosine matrix, all three of 
the other representations are also discussed in detail. 


5.2.1. Direction Cosines 

When a coordinate system B undergoes an arbitrary rotation relative to coordinate 
system A, the basis vectors are related by the equation 

bf = C” A b* (5.2.1-1) 

where the superscripts are coordinate system identifiers, not indices. The matrix of direc- 
tion cosines C DA is orthonormal such that 

C BA C AB _ C AB C BA _ A (5. 2. 1-2) 

It should be noted that the form of the matrix of direction cosines used in this manual is 
the transpose of that developed in reference 21. 

Similarly, with this notation it is easy to show that a basis change for any kinematical 
vector can be performed by changing the subscript and multiplying by the matrix of 
direction cosines for the bases. 

V B = C ba V a (5.2.1-3) 

Note that for kinematical vectors the superscripts are unaffected by these operations. 


5.2.2. Euler Rotations 

If coordinate system B , initially coincident with A , rotates about a unit vector e fixed 
in A by an angle 8 (fig. 5) then the matrix of direction cosines can be written as 

C BA — A cos 8 + e A eA^ '(1 — cos 8) — e A sin# (5. 2 . 2-1) 


where 


_ - u A 

c At — fL * hi 


(5. 2. 2-2) 
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Figure 5. Euler rotation. 

5 . 2 . 3 . Tait-Bryan Orientation Angles 

' B 

Consider two coordinate systems A and B with coincident basis vectors and b { . 
Let the orientation of B with respect to A change as follows (fig. 6 ): 

B 

( 1 ) Perform an Euler rotation of D about e = bj ( j — 1 , 2 , or 3) by an angle By, 

A B 

( 2 ) Perform an Euler rotation of B about e = b k (k — 1 , 2 , or 3, k ^ j ) by an angle 0*.; 

A J5 

(3) Perform an Euler rotation of B about £ = b t {l = 1, 2 , or 3, l ^ k, l ^ j ) by an angle 
01 - 

The final orientation of B relative to A depends both on the magnitudes of 6] , 0 2 , and 0 3 
and the sequence j-k-l. Details of this type of transformation may be found in reference 21 
where Tait-Bryan angles are classified as orientation angles of type body-three. For the 
rotation sequence 1-2-3 the matrix of direction cosines is calculated as follows: 

C 3 S 3 0 C 2 0 — '52 10 0 

C BA -- -33 c 3 0 0 1 0 0 cj 3 j 

0 0 1 s 2 0 c 2 0 3j Cj 

(5. 2.3-1) 

C 2 C 3 Ci3 3 + S]S 2 C 3 Sj3 3 C]3 2 C3 
= ~ C 2 S 3 C 1 C 3 — 5152^3 5 1 c 3 — c 1 5 2^3 

5 2 ^ 1 ^2 c l c 2 


where 


Cj = cos Qi 
— sindi 


( 5 . 2 . 3 - 2 ) 
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Figure 6. Tait-Bryan orientation angles (1-2-3). 


5.2.4. Euler- Rodrigues Parameters 

For two coordinate systems A and B, three parameters (pi = 2 e A i tan(|) may be used 
t.o describe a change in orientation (ref. 21). The values of (pi herein are scaled by a factor 
of 2 relative to the Rodrigues parameters presented in reference 21, so^that for infinitesimal 

values of (pi = (pi, the rotation may be regarded as a vector <^6; = 4>ib { with C BA - A - <p. 
The matrix of direction cosines is then simply 


iBA 


(1 - ffi)A + 


1 + ^ 


(5.2.4-1) 


The angular velocity of B relative .4, expressed in the B basis, can be written as 


n% A = 


(a - 

1 + 


(5. 2. 4-2) 


These relations contain no trigonometric functions and are easily expressed in a shorthand 
matrix notation. Furthermore, a simple inverse transformation exists s ' that given C , 
the values of (j) may be obtained from 

2tijkCfk A 


<Pi 


1 + C BA 


(5. 2. 4-3) 


where C BA is the trace of C BA . Given <p and il% A , '<p can be obtained from 

^ = (A + | + ^)Sl % A (5.2.4 4) 

Transformations between Euler-Rodrigues parameters and direction cosines (or angular 
rates) are very simple relative to the transformations required for Tait-Bryan angles. 
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5.3. Angular Velocity and Virtual Rotation 

The measure numbers for the angular velocity of coordinate system B relative to 
coordinate system A expressed in the A basis, il BA , may be determined from the addition 
theorem discussed in reference 27. They can be related to the time derivative of the matrix 
of direction cosines as follows: 


C BA = 


n% A c BA = -c 


BA n BA 


(5.3 1) 


By virtue of the Kirchhoff kinetic analogy (ref. 28), C BA in Eq. (5.3-1) may be replaced 
with fiC BA , and 0> BA with The expression for the components of virtual rotation 

of D in A then becomes 


sc BA = -S^ B C" A = -C 


BA 


BA 


iBA 


BA 

Ha 


(5.3-2) 


The corresponding virtual rotation vector Sip BA is used in determining the virtual 
work due to applied moments. The components of virtual rotation may be obtained from 
any expression involving the angular velocity in a manner identical to that used to obtain 
equation (5.3-2) from equation (5.3-1). 

Similarly, infinitesimal perturbations of the rotation vector can be obtained by sub- 
stituting C BA for C BA and 0 BA for (l BA in equation (5.3-1). 

C BA = -flg A C BA = -C BA 0 BA (5.3-3) 


5.4, Velocity, Acceleration, and Virtual Displacement 

Velocity and acceleration vectors are obtained by applying the superposition theo- 
rems discussed in reference 27. The calculation of the velocity and acceleration vectors 
is fundamentally nothing more than the differentiation with respect to time of a position 
vector in (i.e., relative to) some coordinate system. It is often necessary to determine the 
time derivative of a vector in coordinate system 27, when the derivative is known only in 
coordinate system A. Given an arbitrary position vector R and its first and second time 
derivatives in A, the first and second time derivatives of R in B may be determined from 
the following expressions. 

b ~R = A ^R + n AB X R 

dt— dt - — 


B 


d 2 

dr 


R - 


B^_ ( B 

dt V dt-J 

d ( a ^r + u ab x r) +n AB x ( A ^-R + n AB xr] 
dt \ dt J \ dt ) 

a ^r+ A ^-n AB x R + 2n AB x A ^-R + n AB x (n AB x r) 

dt 2 ~ dt~ ~ dt~ - — 


(5.4-1) 
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The Kirchhoff kinetic analogy (ref. 28), can also be applied to equation (5.4-1) to 
obtain the virtual displacement vector. Time derivatives in B, B j t { )> are replaced with 
B 6( ); velocity vectors in A, A ^( ), are replaced with virtual displacments in A, A S( ); 
and angular velocity vectors U AB are replaced with virtual rotations Sj/ 30 . 

B SR~ a SR + Si^ AB xR (5.4-2) 
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6. SUBSYSTEMS 


As described in Section 3, the physical structure that is being modeled by GRASP is 
broken down into a hierarchy of substructures. Each of these substructures is represented 
in GRASP as a subsystem. Every subsystem in the model is in turn composed of a set of 
components which may include a frame of reference, a set of nodes, a set of constraints, 
and a set of child subsystems. It is the interrelationship among these components that 
allows the construction of the equations of motion for each subsystem. 

6.1. Frames of Reference 

The position of the frame of reference F for a child subsystem relative to the frame 
of reference 5 for the parent subsystem is defined as R FS , and the orientation (direction 
cosines) of the child subsystem’s frame relative to the parent’s frame is defined as C hs 
(fig. 7). Since Newton’s laws apply only in inertial frames of reference, all equations of 
motions must be written relative to an inertial reference frame. Therefore, it is essential 
to have a method of transforming back to the inertial frame from any subsystem frame in 
the model. If the position and orientation of the parent’s reference frame 5 are defined 
relative to an inertial reference frame 7, the inertial position and orientation of any child’s 
reference frame F can be determined from the parent’s reference frame 5 by applying the 
following equations recursively. 


R fi =R fs ^ R SI 
C FI ^C FS C SI 


( 6 . 1 - 1 ) 


In addition to the position and orientation of any reference frame relative to the inertial 
reference frame, it is necessary to know the inertial motion of every subsystem frame. A 
subsystem reference frame may experience accelerations relative to the inertial frame if it 
or any of its direct ancestors is experiencing translational accelerations or rotation motions. 
Thus, if the inertial motion of the parent’s frame of reference S is known, then the velocity, 
angular velocity, and acceleration of the child’s frame F can be obtained from the following 
equations: 

yFl = yFS + ySI + qSI x gF S 

n FI =n SI + n FS ( 6 . 1 - 2 ) 

a fi =a fs + a si + u SI x ( n SI X R FS ) 
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Figure 7. Frames of reference. 

When expressed in the appropriate bases, these equations become (in matrix form) 

vp =c FS {v s FS + vp + n^RP) 

nP=c FS (n s s f + nf s ) l 6 - 1 - 

Ap 1 =c fs (aP + A s s r + nf J n| J izf s ) 

Note that in the current version of GRASP, it is assumed that V FS — — 0. 



Frames also possess six rigid-body degrees of freedom. Thus, while frames are not 
physically attached to the structure, they may move relative to one another in space. In the 
case of steady-state deformations, these six degrees of freedom include three translations 
along the deformed frame basis vectors and the three Eider- Rodrigues parameters for 
angular displacements. The steady-state displacement vector for frame F is 


R 


F* F 


r? F ' F ; F 
— Rf* « 


( 6 . 1 - 4 ) 


The steady-state frame rotations are expressed in terms of <f> F F , and the direction 
cosines of the deformed frame coordinate axes F' with respect to the undeformed coordi- 
nate axes F are written as Cf- . In matrix notation, the steady-state frame state vector 

( 6 . 1 - 5 ) 



For dynamic perturbations about the steady-state condition, the displacement vector 


R f 


' F' 


= R\ 


:rs f 


( 6 . 1 - 6 ) 


The dynamic perturbations of the frame rotations are expressed in terms of infinitesimal 
rotations 8 F „f , for which the direction cosines (ref. 21) are 


C F"F< 


= A 


~F" F' 

9 F "i 


( 6 . 1 - 7 ) 


In matrix notation, the dynamic perturbation frame state vector is then 



( 6 . 1 - 8 ) 


The virtual displacements for the steady-state and dynamic formulations are simply 
variations of the displacement coordinates 


SR 


F'F 


-SRCfb] 


F' 


F'i 


( 6 . 1 - 9 ) 


F t f 

6R r r =6R r F „ f b, 

and the virtual rotations are variations of the rotational degrees of freedom 




t>± F ' F ^V F ,-b\ 


( 6 . 1 - 10 ) 


6rp 


F F 


" it ' :F 


=6± F , F b ; 
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0.2. Nodes 


Nodes are used by GRASP to model the kinematics of a structure, and their degrees 
of freedom are representative of the physical states of that structure. The position and 
orientation of any node is defined relative to the frame of reference for the subsystem in 
which the node resides. Thus, for a node N in a subsystem with reference frame F, the 
position and orientation of N with respect to F are R NF and C NF , respectively. Two types 
of nodes are currently used in GRASP: structural nodes and air nodes. The kinematics of 
these nodes are described in the following sections. 


6.2.1. Structural Nodes 


A structural node represents a specified material point on a structure. Since the 
material point may have up to six degrees of freedom, the structural node also has six 
degrees of freedom. For the case of steady-state deformations, these six^degrees of freedom 

include three translations along the undeformed nodal basis vectors &,■ and three Euler- 
Rodrigues parameters <j>^' N for angular displacements. The nodal displacement vector for 

node N is then 


R 


N'N 


nN'Nt N 

K Ni 


( 6 . 2.1 1 ) 


The direction cosines of the deformed nodal coordinate axes N' relative to the undeformed 
axes N are expressed as C[j' N . Then, in matrix notation, the nodal state vector is 


q/V' = 



( 6 . 2 . 1 - 2 ) 


Note that the nodal steady-state degrees of freedom are referenced to the undeform.ed nodal 
basis, whereas the frame steady-state degrees of freedom are referenced to the deformed 
frame basis. 


is 


For dynamic perturbations about the steady-state condition, the displacement vector 


R n " n ‘ 


= R 


N"N' 

Ni 



(6.2.1-3) 


The dynamic perturbations of the nodal rotations are expressed in terms of infinitesimal 
rotations 9%'- N ' , for which the direction cosines are 


gN"N' 


-N"N' 

= A — 9 Ni 


(6.2.1-4) 


In matrix notation, the dynamic perturbation nodal state vector is then 

r \ 

= my 


qN 1 


(6.2.1-5) 


Note that the nodal dynamic degrees of freedom are referenced to the undeformed nodal 
basis, whereas the frame dynamic degrees of freedom are referenced to the dynamically 
perturbed frame basis. 
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The virtual displacements fur the steady-state and dynamic formulations are simply 
variations of the displacement coordinates 


6R n ' k =iRm N ii 


6R N " N ‘ =6R%i N 'bi 


( 6 . 2.1 6 ) 


and the virtual rotations are variations of the rotational degrees of freedom 


r / N'N L , N’Nt N 

tt 


( 6 . 2. 1-7) 


r / N" N* cj N"N i 


"/u' : AT 
TV t 


6 . 2 . 2 . Air Nodes 

The generalized coordinates representing the axisymmetric flowfield associated with 
a helicopter rotor are introduced into GRASP by means of the air node. The generalized 
coordinates arc defined relative to an inertial frame of reference /, and determine the 
inertial air velocity at a point Q as 

U' QI = -(Uf + n t + R Q a ? 7,1 + R Q a$ ( 6 . 2.2-1) 

~ A 

where fcj is an inertially fixed unit vector and A is a coordinate system whose origin is 
located at the center of the axisymmetric flowfield. The distance from the center of flow r 
can be calculated from 

r 2 =PO +( R V) (62.2-2) 

For the case of steady-state inflow, U* and 7 ^ represent the uniform inflow velocity 
and the radial velocity gradient, respectively. The other two generalized coordinates have 
no physical meaning under these conditions, and therefore are not used. The air node state 
vector for steady-state inflow is then 

(6. 2. 2-3) 



To model dynamic inflow, generalized coordinates U ^ , 7 ^, and 7 ^ represent the 
collective and two cyclic velocity perturbations. The air node state vector for dynamic 
inflow is then 

fv * 

qA" = < 7ia 

l 7 js 


( 6 . 2. 2-4) 
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7. CONSTRAINTS 


The purpose of a constraint is to create a dependency among generalized coordinates. 
In GRASP, the dependencies among the generalized coordinates are used to eliminate 
dependent generalized coordinates in favor of independent generalized coordinates. In 
the following sections, the general formulation of a primitive constraint will be presented, 
followed by the specific applications in GRASP. Then, the composite constraints that have 
been constructed from the primitive constraints will be discussed. 


7.1. Primitive Constraints 

Consider a set of generalized coordinates that are related to one another through a 
constraint. The constraint relationship g may be written in the special form 


<?e, = 9i (?r 9r w ,)» (* - 


(7.1-1) 


Thus, the generalized coordinates related by the constraint can be partitioned into two 
sets: a set to be eliminated, q e , and a set to be retained, q r ■ Using the constraint relation- 
ship, the set to be eliminated can be obtained directly from the constraint functions which 
depend only on the set to be retained. 

The virtual work for the generalized coordinates associated with the constraint is 


m 


N € 


N r 


^^q ei Qei +' 52 tq n Qr i 

»= 1 i=l 


(7.1-2) 


The sum of the generalized forces Q associated with a generalized coordinate may differ 
from zero for two reasons. First, during the process of seeking an equilibrium solution, 
equilibrium may not always be satisfied. In this case, the sum of the generalized forces is 
residual force that is a measure of the error in the approximate solution. Second, even if 
the complete system is in equilibrium, individual subsystems may not be in equilibrium. 
The generalized forces for these subsystems will be nonzero. 

Taking the variation of equation (7.1-1) 


N r 


i= i qr > 


and substituting Eq. (7.1—3) into Eq. (7.1 2) 


N r 


N e 








(7.1-3) 


(7.1-4) 


This relationship is used by GRASP to incorporate the contributions of the generalized 
forces associated with the eliminated generalized coordinates into the retained generalized 
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forces. During calculation of steady-state residuals, the residuals associated with the elim- 
inated generalized coordinates are transformed and added to the appropriate residuals in 
the parent substructure’s residual vector. 

The treatment of constraints for small perturbations about an equilibrium state is 
a little more involved. For this problem, each generalized coordinate is assumed to be 
the sum of an equilibrium value and an infinitesimal perturbation from that value (he., 
q — q T q ). Equations (7.1-1), (7.1-3) and the generalized forces Q can all be expanded 
in Taylor series about the equilibrium value. Noting that equation ( 7 . 1 - 1 ) is valid when 
q — q, expansion of equation (7.1—1) yields 


N r 


dgi _ 


J dq r . 
i= i H 3 

Expansion of equation (7.1-3) yields 


+ * “ ’ (* - . ,N € ) 


N' / 


d 2 9i 


dq rj " dq r . dq rh 


q r h + 


(7.1 5) 


(7.1-6) 


Expansion of the generalized force, Q , for both eliminated and retained terms yields 


/V e N T 

Qe i "b L € . e . q e . + ^ ^ L e . rj q rj , (i = 1, . . . , N ) 

i=i 

(7.1-7) 

N 9 N r 

Qri —Qn 'b ^ ] L rie . q Ci + L riri q r . , (z = 1, . . . , N ) 

j=i >=i 

where the linear operator, L, contains the terms normally associated with the mass, 
damping, and stiffness matrices, — C — K . Note that the minus signs are present 

in the definition of ju because the generalized force is generally regarded as positive on 
the right-hand side of the dynamical equation, whereas the linear coefficient matrices are 
regarded as positive on the left-hand side. 

GRASP calculat es the M, C, and K matrices for a subsystem by adding the contri- 
butions of each of its children. The rows and columns of the child subsystem’s matrices 
correspond to all of the generalized coordinates of the child. The constraints are used to 
eliminate dependent generalized coordinates, resulting in matrices whose rows and columns 
correspond to only the retained generalized coordinates of the child. The matrices elements 
are then added to the elements of the parent’s matrices that correspond to the child’s inde- 
pendent degrees of freedom. The required transformations can be found using the virtual 
work for the subsystem. 
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An expression for the virtual work from small perturbations about the equilibrium 
state may be obtained by substituting equation (7.1-6) and the eliminated and retained 
subsets from equations (7.1-7) into the virtual work expression in equation (7.1-2). 


AT N T 

i=\ 1 

N r 


N r 


dyi , 


TV* 


N r 


+ E 


9q Vj k=l 


Qr h + ■■■ «~ + E LeiCkQrk + ^ j L ei r h q rk 


fc= 1 


fc=l 


N e 


TV 1 * 


+ ^r< <2r ; + X^ L nej q ej + X] L r irj q rj 


t=l 


i=i 


i — 1 


(7.1-8) 


After discarding terms of second or higher order, the expression for virtual work 
consists of a constant part and two first-order parts in q. The constant part is the same 
as equation (7.1—4), except that it is evaluated for the equilibrium state. 


N ! AT , \ 

1 = 1 J = 1 J 


(7.1-9) 


The first linear portion of the virtual work is the single term 

SS'-fes* 5 -) 

7V r N r 




j—\ k-\ 


(7.1-10) 


The matrix K G represents the geometric stiffness associated with the constraint. During 
assembly of the matrices for the parent substructure, GRASP calculates this geometric 
stiffness and adds it to the stiffness matrix in the parent substructure. This extremely 
important term is often overlooked. For instance, a pendulum, modeled as a rigid-body 
mass constrained to rotate about an offset axis (using a screw constraint) derives all of its 
stiffness from this geometric stiffness term. 

The remainder of the linear terms are 


N T A/' / tV* q— \ 

EE^wE^-.)*.+ 


J— 1 fc = 1 


TV 1 * TV r 


TV e 


dg t 


EE U + E^ 

1 * 3 


j- \ fc=l 


<lr k 


(7.1-11) 
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After substituting equation (7.1-5) into equation (7.1-11) for the eliminated perturbation 
coordinates these terms become 


N r N r 

E5>, 

/ N c 

I L rj rj, + ^ 

j~ 1 k- 1 

V i= 1 


YL re , d91 
’ 9q rk 




N* N e 


dgi 

dqr " jtt dg, i 


L ei e t 


dgi 

dq rh 




(7.1-12) 


The quantity within the parentheses in equation (7.1-12) can be thought of as defining 
a new set of A/, C, and K matrices in terms of the retained and eliminated portions of 
the original matrices. GRASP calculates the new matrices and adds their elements to the 
elements of the parent substructure’s matrices. 

The definition of a constraint follows from the specification of the function g. To obtain 
a solution for a system in equilibrium, the matrix must be known. A perturbation 

solution, however, requires both the matrix and the geometric stiffness matrix K G . In 
the following constraint derivations, matrix will be denoted by 7£. 

7.1.1. Fixed Frame 

The fixed frame constraint describes a rigid connection between two frames of refer- 
nice, F and 5. Regardless of the changes in position and orientation relative to inertial 
space that they may undergo, their position and orientation relative to one another re- 
mains constant. The current (child) frame F will have its degrees of freedom eliminated, 
while the degrees of freedom for the superordinate (parent) frame 5 will be retained. In 
GRASP, this constraint is available through the user interface. 

Steady-State. Consider two frames in their undeformed (5 and F) states and in their 
steady-state (S' and F l ) configurations (fig. 8). The degrees of freedom of F (and F ! ) are 
considered to be dependent, while those of 5 (and S') are independent. The frames are 
assumed to be connected such that 


where 

and 

Thus, 


R s ' s + R SF + R FF ' + R F ' S ' =0 (7.1. 1-1) 

rF', s ' = R FS (7. 1.1-2) 

C S'S C SF C FF' C F'S' _ A (7. 1.1-3) 

R F'F = r F'S' + r s's _ r fs 

R% f = C f ' s '{R% s + Rp') - C F ' S 'C S ' S R FS (7.1.1-4) 

C F'F _ gF'S 1 gS'SgSF 
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Figure 8. Fixed frame constraint. 

Consider the virtual work performed in the F' frame at F' . 

6W = (61 4'. f ) T F£: + (6^% f ) T M% ( 7 . 1 . 1 - 5 ) 

The virtual displacements and rotations at F are related to those at S such that 

8R% F = C F ' S '{8R% S -8C S ‘ S R FS ) 

= C f ' s '(6R% s + W S s' S C s ' s R f s s ) 

= C f ' s '(6R s s : s - R^6i, f,' s ) 

so that the virtual work performed in the F' frame at S' becomes 


( 7 . 1 . 1 - 6 ) 


«W =(6R s s ', s ) T C s ' f 'f£ + 


w§: s ) t (c s ‘ f 'm£: + R$. s Ff) 


FS pF' 


( 7 . 1 . 1 - 7 ) 
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Equations (7. 1.1-7) show the contributions of the force and moment acting at F to the 
force and moment at 5 for the steady-state problem. 

Dynamic. Now consider the two frames in their dynamically perturbed (F u and S" ) 
states (lig. 8). The perturbed position and orientation are related by 



& S"S' +rS ‘F‘ +R F'F" +& F"S" 

(7. 1.1-8) 

where 


^0 

fts 

tl 

lo 

bi to 

(7. 1.1-9) 

and 


c tl c l c / K ,r r ;W r.M rpi! c II 

C s 5 C s h C* f C F s — A 

(7.1.1-10) 

Thus, 

r f ;:, f ' = c f " s "(r s s : s ' + r f s :: s ")- c f " s "c s " s, r f ; s ' 

I pit jtW rit< c* t ^ C ^ C , 1 * 

^rf r t S S S qS r 

Taking the variation of both sides yields 

(7.1.1-11) 


6r f : f ' = c f " s "(sr s s : s ' - sc s " s 'r f : s ') 

qi t qi — S * * 

nfl rrlf of O rpf 

= C F s [£i2f,, 5 +^S” (A-0 S , ),Rf, s ] 

(7.1.1-12) 




vF n S' 


8xp 


s" s' 

S" 


The transformation from the F n frame to the S” frame can then determined in terms 
of 7 Z and Ii G . 

r c F " s " -C k 


n = 


R%- 


o 


iF" S“ 


(7.1.1-13) 


where the columns of TZ are associated with variations of the generalized coordinates 
SRg„ s and 5^5" ? » anti the rows are associated with 6Rp„ F and Sip F ,, F . Then, 


K g 


0 0 
0 FgRg* 


(7.1.1-14) 


where the columns of K G are associated with generalized coordinates R%» s and , 

and the rows are associated with SRg„ s and Sjp^„ s . Equations (7.1.1-13) and (7.1.1- 
14) define the constraint formulation for dynamic perturbations about the steady-state 
configuration. 
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7.1.2. Structural Node Demotion 

The structural node demotion constraint describes a rigid connection between a node 
D in the child subsystem which has a frame of reference F, and a node / in the parent 
subsystem which has a frame S (fig. 9). The degrees of freedom for the dependent node 
D will be eliminated, while the degrees of freedom for the independent node I, the super- 
ordinate frame 5, and the current frame F will be retained. In GRASP, this constraint is 
generated internally, and is not available through the user interface. 






Figure 9. 



Steady- State. The governing equations express the displacement and orientation of 
node D in terms of those of I, 5, and F. The basic equations for the deformations come 

fn ’ m l»' D =R di ' + R 1 ' 1 + R IS ' + R S ' S + R SF + R FF ' + R F ' D 


R 

c D ' D =c d ' i 'c i ' i c is, c s,s c sf c ff 'c f ' d 

In matrix notation the basic equation governing displacement is 

R D'V =c df' c f'f c fs c ss ' [, C S-I^cII'rD’ 1 ' + Rl'l) + Rl/ + RS s 'S} _ 

C DF ' {C F ‘ F R f f s + R f f ', f + RW' ) 


(7. 1.2-1) 


(7. 1.2-2) 
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It is also necessary to take the force and moment at D' and find their contributions to the 
forces and moments at I', F' , and S'. For this, the virtual displacement and rotation of 
D' relative to D is required. The virtual displacement is 


SR 


D'D 

D 


- - C DF 's*‘ F ,‘ C f ' f C fs C ss ' [C s ' r (C ir Rf; r + R f 7 ) + R r s f + R% s } + 


F 1 F 


iF f F n FS r^ss' i ?D' l 


r/> 


>/s' 


>s‘s 1 


*DF' siF' F r'FSr'SS' 


s's , 


C fs C ss S* s , [C s j (c u Rf, 1 + R F1 ) + R* s r + R% s \ + 


/'/ 


C DF 'C F ' F C FS C SS ' [C s ''(Wj M C rr Rfl' r + 6Rj r ) + t>R s s : s ] + 


>s's i 


C DF 6 xJj f , F c F> F C FS R fs C ut "6R % ; 


DF' iT >F’F 
F* 


- c df ’6r f : f + c ds '^i: s + c^sr 1 ; 1 - 
C DF '8t F F F R% F + C DS '6ts. S R$ S + C D/ ^/ V''' 


C7 ds '(^|: s - 

Aft's*}' 1 ) 

and the virtual rotation is 

4 C DS 'W% S - C df '6* f : f (7. 1.2-4) 

The virtual work performed in JP 1 contributes the following terms at I' , S' , and F” 


(7. 1.2-3) 


6W = - (Wl£ F )‘ Fp, + (SR%. s y Fg, + (SR}' 1 )' Ff - (M#' + Rf F F^) 4 


iD 1 




iD' 


/'TV jpD' 


.F'F i 


rO' 


)D'FttD' 


(6*§: s ) (Mg! + Rg s F?) + (S*}' 1 ) (M? 4 - *?''>»') 

(7. 1.2-5) 

Dynamic. For dynamic perturbations about the steady-state configuration, the basic 
kinematic relation is used to determine the matrices 72. and K G . The basic equation 
governing the displ» -ement is 


R 


D"D' 




D"r 


+ R 1 1 4 -R' 1 + R'* 4 -R 

>F' F 


>r'i 


is" 


,s" s' 


+ R 


s's 


rFS _ fiF'F _ ftF" F' _ fiDF" _ rD' D 


(7. 1.2-6) 
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+ R 1 " 1 ' 4 -R 1 1 + R 1 * 4 R 
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F' S' 


R f " f ' - R df " - R d ' d 
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Expressing the position of the perturbed state relative to the steady-state position in the 
D basis (fig. 9), the dependent node displacement may be obtained. 


R D"D> =C DF" C F"F' C F'S' C S'S" C S"I C II' r d''i" + 

C DF" C F"F' C F'S> C S'S" C S"I {R I''I' + R ri )+ 
c df" c f'-f' c f's' c s's- {r I s s" +R sy')- 


,DF" r^F" F' r>F' S' r> s> s " < t>1S" , nS" S' - 
C DF* C F-r- c r-3- R r-'3- + 


,DF" (jF"F‘ C F'S’ rfF'S 
iDF" { t>F" F' , t>DF" 


(Rf, F ' + Rtf )-< D 


The first variation is 
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c DF> 

"c F> 

' F 'c F ' s 'bc s '' 

5 " [C s " I C If, C I,J 

ii 

R?,, 

r + 

c 5 " 1 

r ft 

(Ri 

11 + r i ; i ) + 

f 9 f 1 

Rgn -I- Rs» 

S > 


C°F 

" bC F " F ‘ C F ' s ’ [c s ' 

's" c s"/ c //' c /'y 
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c 5 ' 

s "c s " I (R I I " r 

+ r i ; i ) + c s ' s " 

(R^r 

+ r s s : s ‘) 


( 7 . 1 . 2 - 7 ) 


( 7 . 1 . 2 - 8 ) 


,F’S' 


Similarly, the relationship governing the virtual rotation is 

d"d' . , n" i" , c,l"l’ , clI'I , C././S" , c.i.S" s' 


bip 


--si, D " I " + tr i +6* ri + ty! a " +§r 


bip F ' s - bip L ' l ' - bip 


DF' 


blp 


D'D 


which in the D basis becomes 


H D"D' =zC df" c f"f' c f's' c s's" c s"i 6x(; i I + 

C DF'' c F''F' c F'S' c S'S'' H S‘' iS '_ 


iDF' 


btpp’l'l F> 


( 7 . 1 . 2 - 9 ) 


( 7 . 1 . 2 - 10 ) 
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The virtual work done by a force at D n is then 
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and the virtual work done by a moment at D ,r is 

(6<’ d ') t m£'' =(^r- r ' ) T C is "(A §§: s ')c s ' f '( a + C f ')m£'/ + 
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(7.1.2-11) 


(7.1.2-12) 


The matrices F and K G are 


R = 


C DT _(jDIf(D'l' (JDS' _qDS' rD' S' _£iDF' C DF ' Rp', F ' 

0 C D1 0 C DS ' 0 -C DF ' 


(7.1.2-13) 


where the columns of R are associated with variations of the generalized coordinates 
6R],, 1 , Sift],, 1 , 6Rg,, s , 6xj> g„ s , 6R f „ f , and Sif> F „ F , and the rows are associated with 
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6R and Then, 
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(7.1.2-14) 


where the columns of are associated with 

w, r s s : s ', e s s : s \ and ef, F 
bxpf,', 1 ', 6 r s s ; s ', sr f f ; f \ an d w f : f ' 


perturbations of the generalized coordinates 
\ and the rows are associated with SRp, 1 , 


7.1.3. Screw 

The screw constraint describes two nodes, D and I, that are connected by a mecha- 
nism that permits translation along and rotation about a single axis which is fixed in the 
coordinate system of both nodes. The dependent node D will have its degrees of freedom 
eliminated, while the degrees of freedom of the independent node I wil’ be retained. This 
constraint is available through the GRASP user interface. 

To simplify the derivation, two intermediate nodes located on the screw axis e’ cr will 
be introduced (fig. 10). The “stationary” node 5 is rigidly connected to the independent 
node /, while the “moving” node M is rigidly connected to the dependent node D. Nodes 
M and 5 initially coincide in both position and orientation, but may translate along and 
rotate about the screw axis relative to one another. 
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Figure 10. Screw constraint. 


Steady-State. For the steady-state problem, the equations governing the degrees of 
freedom must be developed, as well as the equation for the contribution of the force and 
moment acting at D' to those at V . The basic displacement and orientation relationships 
for the screw constraint are 


R°' d =R d ' m ' + R M ' 3 ' 4- R 3 ' 1 ' + R 1 ' 1 + R 1S + R sm + R MU 


C D D —C D M M qM’ S' gS'f gfigID 


(7. 1.3-1) 


where M ' indicates a node whose position and orientation relative to S' is the same as 
that of M relative to 5. The position of D ' relative to D in the D basis is then 


Rp D —C DI C 11 ' C rs ' C s '™' C™' M> C M ' D ' Rp, M ' - 

RDM + C DI C II ‘ C 1 ' s< ue’g, r ' 4- (7.1.3-2) 

C^iC 11 ’ R?, 1 ' + R 1 ,' 1 - R 31 ) 

where 

R m ' s ' = ue tcr ' (7. 1.3-3) 


and u is the screw displacement. 


42 



These equations simplify somewhat since C IS = C 1 s — C D M — C DM ^- A. C M * 
can be easily expressed as an Euler rotation, given the screw rotation 9 , and , the screw 
axis unit vector. The virtual displacement is then 


SR 


D'D 


, D =C Dl 6C II 'c I ' s, C s ' lit 'C Iit ' M, C M,D 'Rf j : M ' + 
c D r C C 1 ' s ' C s ' ' S C* 1 ' ' M ' C M ' D ‘ R% M ' + 
C DI SC II 'c 1 ' s ' es c /u + C DI C n 'C I,s ' e‘s c , r ' 6u+ 
C DI (SC IJ ' Rf, 1 + SR 1 / 1 ) 

=c DI [SRj 1 - (. r?' m ' + mr + Rf I ')^I I + 

e‘j cr ' Su + ej er 'R?' M ' 60} 
and the virtual rotation is given by 

S*g D = C Dt (61>}' 1 + e} cr '60) 


(7. 1.3-4) 


(7.1. 3-5) 


The virtual work at the screw connection and at I' due to a force and moment at D 


is 


SW -(SRj T )' Ff' + {SiliftfiR?' 1 ' FP' + MP')+ 
8u(e*r') T FP' + + MP') 


(7. 1.3-6) 


Dynamic. For the dynamics of the constraint, the equations governing the degrees 
of freedom are used to to find the matrices 7 Z and K G . Consider the nodes and the screw 
axis in their perturbed states (fig. 10), an infinitesimal perturbation from their steady-state 
positions and orientations. The basic equations are similar to those of the steady-state case. 


R 


D"D' _ r D"M' 
yl"l 


R 




+ R 


M" S' 


+ R 


s" i 1 


+ 


R 1 1 + R 


I'D 1 


(7. 1.3-7) 


=8V’ M ' +6* M " a " +S+*"*' 


6tl’ S 1 + + 8ip_ 


D" D' c.i.D" M" 
I"S' , c.i.l'D' 


+ 


The first, third, fourth, and sixth terms are zero in both equations. Proceeding as above, 
and noting that 


>,I"V 


=A - C Vl 9 l P r C lV = C v \ A - Of’ r )C 


/'// 




SC 


f,f =-C I ' T {A-6 r 1 " I ')6ri> T i‘ r C 11 ' 


(7. 1.3-8) 


= - Me&'C 
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the 72 and K G matrices are then 


72 = 


c di _ c dirD'i' c DI e , I cr ’ C DI e} cr ' R?' m ' 


iDI 


'Dinner 

e l 


(7. 1.3-9) 


where the columns of 72. are associated with the variations of the generalized coordinates 
f>R[ 1 , (nftj 1 , 8u, and 60, while the rows are associated with SR G D and Sipp D . 


K 


G 


0 0 

0 -FP'R?' 1 ' 

0 — (ej cr ' ) T Fp' 

(Rf' M ‘) T e > , cr ’ Ff' 
— (ej cr ' ) T Mp 


0 0 

— ej cr FP -FP'R?' M 'ei cr ' 

0 0 


0 


( J Rf' M ') 7 eJ cr 'eJ cr > / £ ’' 


(7.1.3-10) 


where the columns of K G are associated with the perturbations of the generalized coor- 
dinates R\ 1 , Op 1 ', u, and 0 and the rows are associated with 6Rj 1 , Syj 1 , 8u, and 
60. 


7.1.4. Copy 

The copy constraint describes the relationship between generalized coordinates that 
are common to both parent and child subsystems, but are otherwise unconstrained. This 
situation most often exists when unconstrained generalized coordinates in the child sub- 
system are passed up to the parent subsystem. This constraint is not available through 
the GRASP user interface. 

Steady-State. From equation (7.1-1), the constraint relationship between the child 
subsystem generalized coordinates q Ci that will be eliminated and the parent subsystem 
generalized coordinates q ri that will be retained can be written as 

ge,=9r 4 , (i = l,.,.,N) (7-1. 4-1) 

Therefore, the calculation of the contributions of these generalized coordinates to the child 
subsystem state vector involves only copying the values of generalized coordinates from 
the parent subsystem state vector into the child subsystem state vector. 

The variation of q ei is then 

6q ei =6q ri , (i = 1, • • • , N) (7.1.4-2) 

When this expression for Sq Gi is substituted into equation (7.1-2), the virtual work is 

N 

*W = £^ r . (Qrj +Qe; ) (7.1.4-3) 

i=i 
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Thus, to assemble the residual vector for the parent subsystem, the contributions of these 
generalized forces from the child subsystem are added to the generalized forces from the 
parent system. 

Dynamic. The derivation of the constraint dynamics follows a similar vein. First, 
the perturbed generalized coordinates, the variations of the generalized coordinates, and 
the generalized forces are expanded in Taylor’s series. 



Qei — Qrj i (* — !»•• 

.,1V) 


(7. 1.4-4) 


he, = 8 Qrj (* = !>• 

..,1V) 


(7. 1.4-5) 


/V N 





~ Qei + ^ei'j Qej + ^ ^ eir / ^ r > ’ 

(* = !,•• 

-,iV) 



i=\ i = i 



(7. 1.4-6) 


N N 




Q n 

-Qn + y Lriej qej + ^ L ri rj Qrj , 

(* = !»•• 

.,1V) 



j= 1 j = 1 





When these expressions are substituted into equation (7.1-2), and the resulting expression 
simplified, the virtual work is written as 


N 


N 


6W = y 8q 


t-1 


Q e . f Qr . + {Laej + L ei r } + L ri ej + q rj 


j= 1 


(7. 1.4-7) 


Since the q ei generalized coordinates exist only in the child system, and the q r< exist only in 
the parent system, I ejr . and L e{ri are null. The 7 l matrix is, therefore, an identity matrix. 
For small perturbations about the steady-state solution, the coefficients in the rows and 
columns associated with the copied generalized coordinates in the child subsystem dynamic 
matrices (M, C, and K ) are simply added to coefficients in the corresponding rows and 
columns of the parent subsystem dynamic matrices. The geometric stiffness matrix K G is 
null. 


7.1.5. Prescribed 

The prescribed constraint is used to describe the permanent deformation of a partic- 
ular generalized coordinate. This constraint is trivial, because the steady-state value is 
constant. In GRASP, the prescribed constraint is available through the user interface for 
nodal degrees of freedom. 

Steady- State. Following the derivation of a general constraint, consider a child 
subsystem that has N e generalized coordinates q ei ,i = 1, . . . ,N e . For this constraint, one 
of those generalized coordinates (e.g., q ei ) has a prescribed, constant value. 


q ei = constant 

q e{ ~9i ( qr 1 5 • • • ) ) ) (* 2, . . . , iV ) 


( 7 . 1 . 5 - 1 ) 
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The total virtual work is 


N c N r 

= ^2 S 9 e iQei ^ SqriQn 

t=I 1=1 

N* N r 


(7. 1.5-2) 


~ ^ ^ ^Q € iQei ^ ^ ^? r *Qrj 
*=2 «=1 

since the variation of the prescribed generalized coordinate 8qe\ is zero. Therefore, this 
generalized coordinate makes no contribution to the virtual work of either the child or par- 
ent subsystem. In practice, degrees of freedom that are prescribed in the child subsystem 
may be eliminated from the parent subsystem state vector. 


Dynamic . The derivation of the dynamic constraint equations for small perturba- 
tions about the steady-state solution proceeds following equations (7.1-5) through (7.1-7). 
The only difference is that in equations (7.1-5) and (7.1-6), i = 2 , When these 
expressions are substituted into equation (7.1-2), 


/V f / 7T Nr ^r~ 


N € 


N r 


i—2 j—\ 


d( b' 3 " 0q rj dq rk 


Qa £ Leie k qe h + ^ L €i r h Qr k 


k - 2 




N r / N* N r 

+ ^ J ^Qr i I Qn + ^ L ri€i q ej + ^ L rirj q rj 

i=i \ i—2 j=i 

(7. 1.5-3) 

From this equation it can be seen that the contributions to the virtual work are the 
same as for the general case, with one exception. The rows of L ei€k and X e . rfc , and the 
columns of L €ieh and J L r . Ci , associated with the prescribed generalized coordinate have 
been eliminated. This is equivalent to removing the appropriate rows and columns from 
the M, <7, and K matrices that are passed up to the parent subsystem. 


7,1.6. Copy Air Mass 

The copy air mass constraint is the constraint used to transform the air mass gener- 
alized coordinates and forces between child and parent subsystems. This constraint is a 
clone of the copy constraint, specialized to copy only the four air mass degrees of freedom. 
Due to the fact that the air mass degrees of freedom are defined in an inertial frame, and 
need never be transiormed out of that frame, the generalized coordinates and forces are 
simply copied. The copy air mass constraint is not available through the GRASP user 
interface. 


7.1.7. Periodic Frame 

The periodic frame constraint describes the relationship between a superordinate (par- 
ent) frame 5 and three or more identical, child frames Fk (for k — 1,2,..., b) rigidly at- 
tached to 5. Frames Fk are located at equally-spaced, azimuthal intervals about an axis 
fixed in 5 (fig. 11). The origin of S is located on the axis of symmetry, while the origins 
of the Fk may be located elsewhere. In GRASP, the periodic frame constraint is available 
through the user interface. 
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The derivation of the periodic frame constraint is very similar to that for the fixed 
frame, except that it is assumed here that there are b identical frames spaced at equal 
azimuthal intervals around an axis. The quantity is independent of k, and C — 

C FlS Tk where 

T fe = T 0 + T c cos (pk + T, sin^fc (7. 1.7-1) 


and where 


T o 


1 0 0 
0 0 0 
0 0 0 


T c = 


0 

0 

0 


0 0 
1 0 
0 1 


T, = 

*k = Y<<=- 


0 0 

0 0 

0 -1 

1), k 


0 

1 

0 

= 1 , 2,. ..,6 


(7. 1.7-2) 


(7. 1.7-3) 


The fixed frame equations can be easily modified to account for this configuration. It will 
be assumed that the axis of symmetry for the periodic, child frames is b t . 
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Steady- State. For the steady-state problem, the equations for the deformed position 
and orientation of any one of the child frames can be written as 


R FLFk =C F ' s 'T k R s s ', s + R F j, S ' - C F '^R F F " k S 


;C F ' s 'T k R% s + R F J: S ' - C ^ 5 ' T k C s ‘ s T F C SF 'R F f 


(7. 1.7-4) 


c nn =C F[s' TkC s's T T C sF l 


In order to make the left-hand sides independent of fc, let Rg,f = — 0 and S ~ 

<f > 3 s — 0 . Since the right-hand sides are equivalent for all fc, all T* can be set to To to 
simplify the equations. The virtual displacements are then 


SRp! Fk =C f ' s 't 0 6R s s ', s - C F ' S ' T 0 6C s ' s T F C SFl R F 1 s 


(7. 1.7-5) 


5'S 


-C F J s 'To6iZ|! s + C f ' s 'T 0 8iI> si C s ' s T F C SFl R F 1 s 

and the virtual rotations are 

64>2 Fk - C F ' S ' To6i/>g. s 


(7. 1.7-6) 


The virtual work at 5' due to the b sets of forces and moments acting at F £ is therefore 
6 


[(**$*) f^+{w f f i Fk ) m f f 


fc — 1 


-b{(6R s s : s ) T T F C s '^F F i 4 - (6^i: s ) T [TjC s ' F '^M F i- 

( tJ c s ' f * f f { )~c s ' s t f c SFi s ] } 


(7.1. 7-7) 


Dynamic. For small perturbations about the steady-state solution, the perturbed 
position and orientation (fig. 11 ) of any one of the child frames F k is 


R n'n 

u F ,, 


R F i’, S " +C F " s "T k R s s : s ' - C F " s " T k C s " s ' C s ' F ' R Fk , s 


(7. 1.7-8) 


C FL'FL = c F[‘ s" TkC S" S' T T c s' f[ 
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To first order in the perturbation quantities, the virtual displacements and rotations are 


Fj!K =C K 's" Tk \ 6R S”'S‘ + H , s sn s (A _ r s ,r )T?C s ' F *R F F l 




6R\ 




S'tjj 



F['S" 


w 


5"S' 

5" 


(7. 1.7-9) 


~S"5' 

Note that the geometric stiffness matrix will come from the 9 S „ term. 
After substituting into the expression for the virtual work, the matrix 77 is 


n = 


C F " s "T k 

0 


-{T F c s ' F iR Ks ' F '>y 


C F " 5 T k 


(7.1.7-10) 


where the columns of 77 correspond to bR s s ',', s ' and and the rows correspond to 

bR F 'i F ' h and Sip F ', F ' r . The matrix K G is 

t k * 


K g = 


0 


t=i 


(7.1.7-11) 


where the columns of K G correspond to R$t f S and 9$ t , s , and the rows correspond to 
6R F i Fk and 6V* For evaluation of the lower-right submatrix, it should be noted that 

F h r 


£ T?( )T„ = b\T,T( )T. + *!?( )Tc + I TT( )T.\ (7.1.7-12) 

1 

when the expression enclosed in parentheses is independent of k. 

7.1,8. Periodic Node Demotion 

Just as the periodic frame constraint is very similar in concept to the fixed frame con- 
straint, the periodic node demotion constraint has a similar relationship to the structural 
node demotion constraint. In this case, a node belonging to a parent subsystem is repli- 
cated in the child subsystem at b equally-spaced azimuthal intervals about an axis that 
is fixed in the parent subsystem. The periodic node demotion constraint is not available 
through the GRASP user interface. 

The degrees of freedom of the 6 child subsystem nodes D k are expressed in terms of 
the degrees of freedom of TV, S', and I. To visualize this constraint, consider figure 11 and 
imagine a node I associated with frame S and a node D k associated with each frame F k , 
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as in figure 9. The virtual work done at all nodes D k for k = 1,2, ... ,6 is determined at 

F k . The total virtual work is summed for k — 1,2, ... ,6 and determined at 5 and I. For 

. 

this constraint, it is assumed that the axis of symmetry is b x , and that b x also lies along 
that axis. 

Steady- State. The governing equations for the periodic node constraint are derived 
in a. similar manner to those of the structural node demotion constraint. First, let 

C FhS = C FlS T k = C F ' s *T k 

(7. 1.8-1) 

C DhI =C D '" J ' = C DlI T k = C D ' V T k 
where T k is defined in equation (7. 1.7-1) and (7. 1.7-2). When 

Rf h h S =Rp h , S = constant 

(7. 1.8-2) 

R% h J = R = constant 

the positions and orientations of the frames and nodes may be written as 


K u k 


zC D >nc F '> F *C F ' s T k C ss ' [C S / (C 7/ C I Dl Rp? 

Ri r ) + R'/ + R% s }~ 

C DhF '»{C F '>' Fh Rp1 s + Rpi Fh Rp^ F! ‘) 


(7. 1.8-3) 


qD^Di. TkC 1 1 C IS C s s T k C SFx C FkF>h C F ' uDk 

To make equation (7. 1.8-3) independent of fc, let Rg? — 0, C IS = A and let T k — To (all 
displacements and rotations for 5 and I take place along or about the axis of symmetry). 

The virtual displacements and rotations are required in order to calculate the virtual 
work of a force and moment at D r k for all k. 


- = - C^ F ^^ F ^ FiF ‘C Fl5 r u C 5S '[C 5 ' / (C // 'T 0 T C iD 'i2^ / ' 

R t , I ’) + n t / + R i: s ] + 

C D ' F lc F '* Fh C F ' s TvC ss, fy S s, S [C s 'i{ C ii'Tjci D 'R%*\ 

R l /) + R , /+RS s ] + 

c^ncn^c^T.c 55 ' [c s ' I {U I i I c II 'Tjc I ' D ''R D jf+ 

sri'^+sr i : s j - 

C D ^{~H F ;^C F ^R F f + 6Rp, Fk ) 
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(7. 1.8-5) 


=C D ^C F '^C F ^ S T,C SS \C S,I 6^; 1 

c DhF ^Si^ Fh 

The virtual work due to the virtual displacements of each of the nodes D' k is 

tfW = - (6R% Fk ) T F% - + *#*<)+ 

(6R% s fc s,s T 0 T C SF 'Fg k + 

WI: s ) T {C s ' s T f C sf 'M% + [C s ‘ s T F C SF 'C F ^{R D F f k 4- C F * F >R F f+ 

C f > f *C f ' s ToC ss 'rI\ s )Yc s ' s T f C sf 'F°*} + 

{bR I i I ) T C IS 'c s ’ s T o T C SFl F Fh h + 

{f>Tl> I i I ) T (C IS 'C ss T F C SFl Mp* + {cWlfc'W- [R F j; Fh + 
c n F >R F f + c F * Ffc c Fi5 r 0 c ss '(f7|', s + #3?)+ 
cFiFucFiSTC^C^^l'^yc^'c^^C^F^) 

(7. 1.8-6) 

The summation of terms involving the virtual displacements and rotations at 5 and I 
involves only a multiplication by b. The corresponding terms at F k need not be summed 
since only one system contributes to the virtual work there. 

Dynamic. As in the case of the equations for the steady-state periodic node demotion 
constraint, the derivation of the equations for the dynamics is similar to the derivation of 
the structural node demotion constraint dynamic equations. In a manner similar to that 
for the static equations above, let R{^ = 0 and C IS = A. Also, let 


c ns ^cFis' = c F' h ’s" = c FlS T k = C F ' s 'T k = C F " S ' T k 


constant 


c D k I =c D' k r' _ c d';i" _ c DlI T k = C D ' f Tu = C D " T "T k = constant 


(7. 1.8-7) 


nA 

il n 


F k s =nT =*% 


F’J S' 


constant 
= R n ~ constant 


51 


The resulting equation for the matrix 7 Z is 


1Z = 


C D - ! -R^t'c 0 ' 1 C D - S ' -c d - s 'r%- s ' 


Du 

C D -* 


C D - S ' 


c D >n c D -nR j* 


D’F 


c D '* F > 


(7. 1.8-8) 


where 


C DuI =C D - F -C F ' s 'T k C SI 


C D - S =C D - F '-C F 'i s 'T k 


(7. 1.8-9) 


R%- S ' ^'^Ttcno-R* 1 ' 4- C s,i R Fi + R 1 / 

The columns of TZ correspond to bR\ 1 , bt^j F , bRg',' s ' , brpg',', s> , 6Rp-, Fh , and bip^-, Fk , 

D n D* D n D 1 

respectively; while the rows are associated with h an d h ' The coefficients of 

the geometric stiffness matrix K G are then 


6Rj 1 row: 


fit r 

Rj column: 0 


jit ,1 

6j column: 
column: 
6$„ S column: 
R^ h n Fk column: 

F k 

n F hK i 
p rn column. 


0 


0 


b(T F F^y 


0 

T k C IFl FpS 


(7.1.8-10a) 
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rlt t 

Sip j row: 


R\ 1 column: 0 


column: - VT* (ff‘ fl? 1 ' )7i - ^Tj (F?' fl?' 1 ' )T C 

h -Tj(F?'R°'‘)T, 


R%i S column: 0 


(7.1.8-10b) 


i r 1 - n f 


column: - F?') T, - F?’ )T C 


b -Tj(Rp , 'Ff’’')T, - b(TjM f')' 


column: 0 

F k 


O^f" column: TZ R^ 1 ' C if 'F$ + T k C IF '‘M°P 

** it 11 


SR%„ row: 


ilj 1 column: 0 


1 column: 0 


iZf.i 5 column: 0 


0f„ 5 column: —b{TjFf l ) 


Rp k „ Fk column: 0 


(7.1.8-10c) 


*5^ column: T k C IF '^FpP 


K' 
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iZc,, 5 column: 0 

6 (7.1.8-10e) 

0s ft s column: 0 

R^t Fh column: 0 

0^ h ,, Fk column: 0 

F k 
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row: 


rtl rt 

R, column: 


- F^C F ^T k 


0j 1 column: Fp, 1 C FlT R.f l1 Tk 

il|„ s column: — FpSC FlI T k 

0f„ 5 column: Fp, l C Fl1 R® lt Tk 


(7.1.8-10f) 


R 1 !), Fk column: 0 


0pn Fh column: 


- n' - n* F* 

F u i n lf i 1 i 

pt itrf 
r i ' i 


7.1.9. Periodic Generalized Coordinate 

A rotationally isotropic structure consists of three or more identical substructures 
that are spaced around an axis of symmetry at equal azimuthal intervals. The periodic 
generalized coordinate constraint exists in order to transform generalized coordinates that 
belong to the rotationally isotropic structure into the generalized coordinates for a generic 
member of that structure. Additionally, it must transform generalized forces for a generic 
substructure into the generalized forces for the complete structure. In one sense, it is 
simply an extension of the copy constraint for periodic structures. In GRASP, the periodic 
generalized coordinate constraint is not available through the user interface. 

Steady- State. The set of independent generalized coordinates for a rotationally 
isotropic structure may be grouped as collective (q 0 ), cosine ( q c ), and sine (q 3 ) components. 
The generalized coordinates for the fcth generic substructure q k may be written as 

q k = q 0 + q c cos (f>k + q» sm<j>k (7. 1.9-1) 


where <f> k = — !)• The variation of these coordinates is 


8q k = 8qo + 8q c cos (f>k + 8q t sin (f>k 


(7. 1.9-2) 


Given generalized forces Q k , the total virtual work from all of the generic substructures 

8W = = 8q FS ^Qk + 8q c T ^ Q k cos (f>k + 8q, T ^Qk sin <f>k (7. 1.9-3) 

fc=i fc= i *=i fc=1 
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Since the generalized force of a generic structure is independent of k, Q v = Q k and 

6W = Sq 0 T bQ 0 (7. 1.9-4) 

Dynamic , The dynamic perturbations of the generalized coordinates are related in 
the same manner as the variations of the steady-state generalized coordinates. 

qk = qo + q c cos 4>k + q 9 sin((} k (7. 1.9-5) 

Like the generalized forces, the substructure coefficient matrices M, C, and K are inde- 
pendent of k and 

OV = 6q k T {Mj k + C'q k + Kq k ) (7.1.9-6) 

The contribution to the virtual work in terms of the independent generalized coordinates 
is 


n 

SW — ^ ^gu T + Sq c T cos 4> k 4 - Sq, r sin^fc 


Qo 


k- 1 


[M] { q c cos </>*•} + 

q . si« 4>k 


9o 1 I 9o 

[C] { q c cos <f) k > + [K] < q c cos <f> k 
q, sin 4> k J [ q. sin (f> k ) 


(7. 1.9-7) 


b6q 0 ( Mq 0 -f Cq 0 + if go) + ~^gc (Mg c + Cg c + Kq c ) 


+ - Sq, 1 (Mg, + Cg, + Kq,) 


The matrices for the rotationally isotropic structure therefore have three rows and columns 
for every row and column in the generic substructure, and are of block diagonal structure. 


7.1.10. Periodic Air Mass 

The periodic air node constraint describes the transformation of the air node gener- 
alized coordinates and forces between subsystems associated with periodic structures and 
subsystems associated with generic substructures. Since the air node generalized coor- 
dinates describe an induced airflow velocity field that is already axially symmetric, the 
periodicity of the structure has no effect on them. In fact, it is assumed a priori that the 
flow field is interacting with a rotating, periodic structure. This constraint is not available 
from the GRASP user interface. 

Steady-State. When a subsystem is periodic (in the sense that it consists of three or 
more generic, periodic members such as those described under the periodic node demotion 
constraint), the steady-state air node generalized coordinates f7f 4 and 7 ^ are simply copied 
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from the parent subsystem to the child subsystem in a manner similar to the copy air 
mass constraint. During the assembly of the generalized forces, the air node generalized 
forces from a single, generic substructure are simply multiplied by b and added to the 
corresponding generalized coordinates of the parent subsystem. 

Dynamic. For perturbed motions, let the generalized coordinates for the fcth subsys- 
tem be v 

<jk = Tkq', bq k = Tkbq (7.1.10 ) 


where is as given in equation (7.1.7 1) and where 


and 



f fp '\ } 

r spf ) 

HU 

Or 

II 

On 

1 HU J 

y ^^13 J 


(7.1.10-2) 


(7.1.10-3) 


Note that q k and q will not appear in the dynamical equations. The equations now trans 
form in exactly the same manner as the ones in the copy air mass constraint. 


7.1.11. Rotating Frame 

The rotating frame constraint describes a constraint that is very much similar to the 
fixed frame constraint, except that frame F is rotating at a constant angular speed relative 

to frame 5 (fig. 12). The axis of rotation passes through the origin of F and along fc. 
No time-dependent terms are retained in the equations. In GRASP, this constraint is 
available through the user interface. 

Steady-State. In moving to its steady-state, equilibrium position, the axis of rotation 
follows &f. The position vectors R SF and R S F are constant in the 5 and S' bases, 
respectively. The change in orientation is then 

C F ' s \t) = T{t)C F ' s '{0) = C FS (t) = T(t)C FS (G) (7.1.11-1) 


where 

T(t) = To + Tc cos («0 + T. sin(nt) 


(7.1.11-2) 


where To, T c , and T t are given in equations (7.1.7 2). 

The kinematics for the rotating frame constraint are based on the following equations. 


pF'F 

pi 


TC f ' s '( 0 ) (R 


f' s' 

S’ 


+ Rg, s - C s ' s R fs ) 


C F ' F =TC f ' s ’(0)C s ' s C sf (0)T t 


(7.1.11-3) 
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The time-dependent terms in these equations vanish if all of the displacements and rota- 

~F' 

tions of S' relative to 5 are along and about fe 1 (the axis of rotation). Therefore, let 


=R^ = 0 


/F'F _j.F'F _ n 
02 03 — 0 


(7.1.11-4) 


The virtual displacements and rotations of the F f frame are then 




6R£ f = { o } = T„C F ' s '(0) (SR s s '. s + ii,s' S C S ' S Ks S ) 
0 


SVf'' F 


Hf'i 

o } = ToC f ‘ s '{0)6iP§', s 
0 


(7.1.11-5) 
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where the use of only the Ty component of T eliminates the time-dependent terms. 

The virtual work at 5 associated with the force and moment at F yields the following 

T 


contribution at S: 


6W =(tRs’ S ) C SF (Q)T 0 r F£, - 

(Srp s s : s ) T R^ s C s ' F \Q)T F F F : + 


(7.1.11-6) 


W I-' 5 ) c s ' f '(o)t f m f ; 


T ,/F' 


Dynamic. The position and orientation of the perturbed frame relative to the steady- 
state position and orientation are related as follows: 


r f ; f ' =tc f " s "( o) (r f :: s " + r s 5 : s ' - c s " s 'r f : s ') 

C F"F’ =T c f " s "(0) c s " s 'c s ' f \ 0) T t 


(7.1.11-7) 


From these equations, the virtual displacement and rotation may be obtained. To first 
order in the perturbation quantities, 


sr f " fi =rc F " s ''(o) [6r s s " s ' + *4* (a - ~C. S )^f' 5 '] 


(7.1.11-8) 


S i i f ’ f 

o 


R is 


=Tc F " s \o) 6 i> s s ;:' 

where contributions due to geometric stiffness come from the 0 S „ term. The matrix 


-5" S' 


R 


T{t)C F " s "{0 ) -T(t)C F " s "{0)R F ',‘ 

o r(t)c F '' s "(o) 


(7.1.11-9) 


where the columns of R correspond to R^',', s ' and 0f" 5 ', and the rows correspond to 
SR f ,, f and f>ip F ’,, F . 

Since R depends on t, the time- dependent terms must be removed from the final 
transformed equations. This is easily accomplished by taking the time-averaged value 
of the transformed equations. The oidy contributing (i.e., nonzero) terms then are the 
constant terms, the cos 2 (flT) terms, and the sin 2 (flT) terms. In addition, since R depends 
on t, terms from matrix M will contribute to M, C, and K in the transformed equations 
and C will contribute to C and K by virtue of the following relations 


(7.1.11-10) 
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(7.1.11-11) 


F n F* 

Rptt 

• F" F* 


= R 


T>S"S' 
K S U 


nS"S' 

V S " 


+ 2 R 


. s" s' 

R S " 

•5" S' 


r pit / ' ~ " V^5/ 

K G i the geometric stiffness matrix, is 


K u = 


+ R 


s n s' 

R S n 

’ s n s 1 


0 0 
0 lC s " F \Q)TTF£;}~Rr; s ’ 


(7.1.11-12) 


where the columns of K G correspond to and , and the rows correspond to 

SRg* and fyg*. 


7.1.12. Rotating Node Demotion 

The rotating node demotion constraint describes the relationship between two nodes, 
one of which is located in a rotating frame of reference and the other in a nonrotating frame 
of reference (llg. 13). This constraint combines many of the characteristics of the rotating 
frame and structural node demotion constraints. It is assumed that the child frame F is 
rotating about a fixed axis at an angular speed ft, and that a dependent node D is defined 
relative to that rotating frame. The parent frame 5 is stationary (relative to F) and an 
independent node I is defined relative to 5. The rotating node demotion constraint is not 
available through the GRASP user interface. 

Steady- State. The governing equations for the steady-state condition are similar to 
those for structural node demotion except that 

C FS {t) =T{t)C FS { 0) = C F ' s \t ) = T(/)C f ' 5 '(0) = C DI (t) 

=T{t)C DI {0 ) = C D ' r '(t) = T{t)C D,I \ 0) 

In addition, R$ S and R®, 1 are constants, and 

T(^) — To 4" T c cos(m) 4" Tj sin(fl/) 


(7.1.12-1) 


(7.1.12-2) 


The governing equations describing the deformed position and orientation of the de- 
pendent node are then 


rD'd ^ c df' c f'f tc fs [ 0 )c ss' [ c s>j {c ii' r d;i' + £/'/) + R i? + R S S S} 


C DF' [c f ' f TC fs (0)R% s 4- R F ' F + Rp' F } 


F'F 


,DF' 


(7.1.12-3) 


C D ' D =tc d ' v 


(0)C F 'C /5 'C S ' S C SF (0)T T C FF 'C F ' D 


In order to be independent of t, let T — To and choose Rp, F — 0 and C DF — A. Thus, 
only displacements along and rotations of D about the axis of rotation can be nonzero. 
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Figure 13. Rotating node demotion constraint. 


The virtual work at D' due to a force and a moment acting at D' is determined in terms 
of the virtual displacement and rotation 


6R°' d = - C df 'W F f , ' F C f ’ f T 0 C fs (Q)C ss ' +C DF 'C F ' F TC FS (0)C ss 'R° s 'wi S + 

C df ' C f 'f ToC f-\o) C ss ' [C s ' I {S^ r i I C II 'R?' 1 ' + SR 1 / 1 ) + 6R S S ', S ]- 

C DF ' [ _ ty F F F C F ' F T 0 C FS {0)R% s + 6R f ', f ] 

=C DF T 0 C FS (0)C SI 6i> I I ' 1 + C DF r 0 C F5 (0)C S5 '6V»Ir' s - C df 'H'f F 


(7.1.12-4) 
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The virtual work done at D* is then 

SW = - (6R f f : f ) T F°: - (Si,£ F ) T (Mg + R# f F#)+ 

(6r§ , , s ) t c s ' s c sf (o)t 0 f£ >, + 

(6ij>§', s ) T [C s ' s C SF (0)T?Mi?' + R%,' s C s ' s C sf (0)T f f£ > '}+ 


(7.1.12-5) 


(SR 1 / 1 ) C IS C SF (0)T {i Ff > ' + 

(Si)’ 1 / 1 ) T [C 1 s C F s (Q)T f M ft' 4- Rj > i C 1s C sf (0)T f F f ] 


Dynamic. The governing equations are similar to those for structural node demotion 
except that, as in the static case, R^F — 0 and C DF = A. The governing equation for 
the position of the perturbed dependent node is 
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Solving for the dependent node displacement, 
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D " D ’ -nD"V +R l" V +R I- 1 +R IS +R 
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- R df - R l 


S H S' 
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By referring the displacements to the D basis this vector relation becomes 

1°' 

c o 


R 0"D' _ C DF" C F" F' gF' N' £N' S' C S‘ S" gS" I gll' gl' I" rD'* I u + 
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q DF" gF"F' qF' N' qN' S' R N' S' 
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DF" (jF"F' r^F'N' pF'/V' 


pFWV' r>DF" ( r>F n F' , p U£ ”\ ryU 

~r 


DF' 


jD'D 


The virtual displacement is then 

8R%' d ' =C DF ”C F " F 'C F ' N 'C N ' S 'C ^"c^^Rf'^ 

C 0F" C F"F' c^^C^^C^^SRiy - C df " 6R f : f ' + 

C DF " C F " F ' C F ' N ' C N ’ s ’ C s ' s " C s " 1 C 11 ' 6C 1 ' r " Rfn T " + 
C UF ” C F " F ' C F * N ' C N ' s ' 6C S ” s ' [C s ' I C lF C jF ’Rp” I " + 

c s " i (R i I ,,if + R 1 ; 1 ) + R j s r + r s 5 : s '}+ 

C DF "sc F ” F 'c F,N, [C N ' s 'c s,s ''C sni C II 'c I ' r "Rf" I '' + 

c N's'c sts "c s " I {Rf' 1 ' +Rf*)+ 


iN's' r^s's" < nis" i p 5 " s' \ nN’s' r n's‘ _ 


(Rs* FR^n-c 1 


F'N’ 

N' 


(7.1.12-6) 

(7.1.12-7) 

(7.1.12 8) 


(7.1.12-9) 


(7.1.12-10) 
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and the virtual work done by the force at D' is 
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(7.1.12-11) 


The virtual rotations can be easily obtained from the variation of the direction cosine 
relation 
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Upon removal of the tilde the virtual rotations are 
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(7.1.12-13) 


(7.1.12-14) 
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and the virtual work done by the moment is 
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(7.1.12-15) 

since C F ^ — T, let $ = C F F , as in the structural node demotion constraint. The 
time-dependent terms can be eliminated when R F , F = 0 and C DF = A. 

Combining the virtual work due to the force and the moment at D" , the matrix 77 
can then be calculated. 


77 


c df' tc f'i _c df ' TC f' i rD' /' c df' tc f's' 
0 C df 'TC f ' 1 * 0 


(7.1.12-16) 


-c df ‘tc f ' s ’r?' s ' 


l S' 


iDF' jp^F' S' 


C DF ' 0 

0 - c DF ' 


where the columns of 77 are associated with 8R F i ' , Sijjj" 1 ' , 6Rg,', s ' , Sipgn S ' , 6R F " F , and 
t^ F ,, F , respectively; and the rows correspond to 6Rjj' D ' and Sip® " D ' . The coefficients of 
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the matrix K G are 
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i" i' 
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row: 



row: 


Rj 1 column: 0 

0j 1 column: 0 


R s s „ 5 column: 
#f" s column: 
Rp ,, f column: 
9 F ", F ' column: 

rft rf 

R, column: 

rll r> 

9f column: 

(jf I c/ 

R$n column: 

rll qf 

Ogt, column: 
Rp„ F column: 
6 F ,, F column: 


0 (7.1.12-17a) 
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row: 
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row: 


r/f r/ 

Rj column: 0 

0j 1 column: 0 
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Rs'" S ' column: 0 (7.1.12-l7e) 
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(7.1.12-17f) 


it ’ ^ ^ pr ^ 

6 f ,, column: 0 

7.1.13. Rotating Generalized Coordinates 

The rotating generalized coordinate constraint relates generalized -oordinates in one 
subsystem to the corresponding generalized coordinates in another subsystem that is ro- 
tating at constant angular speed relative to the first. This constraint is often applied to 
subsystems that contain periodic structures. This constraint is not available through the 
GRASP user interface. 

Steady-State. The general form of the transformation from rotating to nonrotating 
coordinates is 

, fi = r, N (7.1.13-1) 

where 

T = T 0 f T c cos fit + T, sin fit (7.1.13-2) 
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and q n is a set of generalized coordinates in the rotating subsystem that corresponds to a 
set of qn generalized coordinates in the nonrotating subsystem. 

In order to make the transformation equations independent of time, let T = T 0 and 
ftq R T — TySqw. This eliminates any generalized forces of the lateral (cosine or sine) type. 
Then, the virtual work is 

6W = qlQ R = qJ r T?Q R (7.1.13-3) 


Dynamic. In the rotating system, the virtual work can be written as 


<5W = q R {M R q R + C R q R + K R q R ) (7.1.13-4) 


where 

qR =Tq N 

q R =TqN + Tq N 


(7.1.13-5) 


—Tq^ + 2 Tq N + Tq N 

Substituting these relations into equation (7.1.13-4), the virtual work can be obtained in 
terms of the generalized coordinates of the nonrotating system. 

m - qJfT T \m R T'4 n + ( CrT + 2 M R T)'q N + {K r T + C n T + M R T)q N ] (7.1.13-6) 

Thus, the C N coefficient matrix (in the nonrotating system) depends on C R and M R , and 
K /v depends on K Rl C R , and Mr. 


7.1.14. Rotating Air Mass 

The rotating air mass constraint transforms the air node generalized coordinates and 
their associated generalized forces between a rotating subsystem and a nonrotating sub- 
system. As in the other air mass constraints, the air node generalized coordinates are not 
transformed out of the inertial frame of reference. The rotating air mass constraint is not 
available throught the user interface in GRASP. 

Steady-State, oince only U* and 7 j* r are valid coordinates in the steady-state prob- 
lem, and both are rotationally symmetric, they are treated in exactly the same manner as 
in the copy air mass constraint. 

Dynamic. For a set of dynamically perturbed air node generalized coordinates, let 
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(7.1.14-1) 
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and 


(7.1.14-2) 


Then 
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qR =Tqw 
qR = Tqw 
qR =Tq N + Tq N 


} 

Win 

^ 13 N } 


(7.1.14-3) 


The virtual work for the rotating subsystem is 

6W = ql{M'l R + C'q R ) (7.1.14-4) 

which, when transformed into the nonrotating subsystem, becomes 

iW = q lT r [(A/T<j w + (CT + MT)ij N ] (7.1.14-5) 

The C coefficient matrix for the transformed (nonrotating) subsystem therefore depends 
on the M and C coefficient matrices from the original (rotating) subsystem. 


7.2. Composite Constraints 

In general, a composite constraint is a constraint that is built, np out of one or more of 
the primitive constraints that have been described in the previous sections. The bundling 
of primitive constraints into a single constraint is primarily done for the convenience of the 
user There are many times that sets of constraints must be used together, and it makes 
sense to combine them internally. In the following sections, the composite constraints that 
have been constructed from the set of primitive constraints in GRASP will be descri let - 
All of the composite constraints are available from the GRASP user interface. 

7.2.1. Aeroelastic Beam Connectivity 

The purpose of the aeroelastic beam connectivity constraint is to provide a means for 
attaching an aeroelastic beam element to a GRASP model. The element subsystem for 
the aeroelastic beam consists of a frame of reference, a root node, a tip node, and an air 
node, all of which must be connected to their counterparts in the existing portion of the 
model. Therefore, the aeroelastic beam connectivity constraint must contain a fixed-frame 
constraint (for the frame), two structural node demotion constraints (for the root and tip 
nodes), and a copy air mass constraint (for the air node). 

In the definition of the aeroelastic beam connectivity constraint, the position and 
orientation of the dependent, element root node R relative to an existing, independent node 
I n must be provided. The position and orientation of Ir relative to its subsystem frame of 
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reference (the superordinate frame Sr) is known from the definition of Ir. Therefore, the 
position and orientation of the dependent, root node relative to the superordinate frame 
can be calculated. 

(jRSr .— (jRIr(jIrSr 


R r ^ r =c Snln R^ R + 


Jr 


l s r 


( 7 . 2 . 1 - 1 ) 


After locating the parent subsystem of the element subsystem in the system organization 
tree, the position and orientation of the parent frame relative to the superordinate frame 
can be calculated. Since the element frame and the element root node are coincident, 
the position and orientation of the element frame relative to its parent frame can then be 
determined. 

C FP —(jRP _ (^RS r( jS r P 


( 7 . 2 . 1 - 2 ) 

Rp P =R$ P = C PSr (R r ^ r + Rs* P ) 

With this information, the fixed-frame constraint can be defined. In addition, all of the 
position and orientation information is available to define the structural node demotion 
constraint for the element root node. In those cases where the superordinate frame is 
not the same as the parent frame, it is necessary to create copies of the independent and 
element root nodes in each of the subsystems leading to their nearest common ancestor. 
These nodes are chained together by a series of structural node demotion constraints. 

The position of the element tip node T relative to the root node is defined as R tr = 

£b 3 , and the orientation C' 1R is defined as an Euler rotation of magnitude 0 f £ about 6 3 . 
After the position and orientation of the root node relative to the independent tip node 
It has been calculated, the offset of the element tip node from the independent tip node 
can be determined. 


R]* T = C 1tSt (C StR R^ r + R RSt - iZ' rSr ) 


qTIt ^(jTRfjRSr (jStIt 


( 7 . 2 . 1 - 3 ) 


At this point, the structural node demotion constraint for the element tip node can be 
defined. In those cases where the superordinate frame is not the same as the parent 
frame, it is necessary to create copies of the independent and element tip nodes in each of 
the subsystems leading to their nearest common ancestor. This creates another chain of 
structural nodes, all connected together by structural node demotion constraints. 

If the beam element is to be connected to an air mass element, the position and 
orientation of the element subsystem relative to the corresponding air node is calculated. 
Then, the copy air mass constraint is defined. In those cases where the air node is not 
defined in the parent frame, it is necessary to create copies of the independent and element 
air nodes in each of the subsystems leading to their nearest common ancestor. This creates 
a chain of air nodes, all connected together by copy air mass constraints. If, however, the 
beam element is not to be connected to an air mass element, the four nodal air mass 
degrees of freedom are constrained out using prescribed constraints. 
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7.2.2. Air Mass Connectivity 

The purpose of the air mass connectivity constraint is to provide a means for attaching 
an air mass element to a GRASP model. The air mass element subsystem is unusual in 
that the frame serves only to establish the position and orientation of the element relative 
to the remainder of the model. Therefore, while the frame does exist and does need a 
frame constraint, it has no frame degrees of freedom. The air mass connectivity constraint 
is then made up of a fixed-frame constraint, a copy air mass constraint, and one or more 
prescribed constraints. 

In the air mass connectivity constraint it is assumed that the independent air node /, 
the dependent (element) air node A, and the element frame F are all coincident. 


C tA =C AF = A 

rIA = RAF = 0 


(7.2.2-1) 


After locating the parent of the element subsystem, the position and orientation of the 
parent frame relative to the superordinate frame can be calculated. Since the position and 
orientation of the independent air node relative to its subsystem frame (the superordinate 
frame 5) is known, the position and orientation of the element frame relative to the parent 
frame can also be calculated. 


C FP =C IP = C I5 C SP 

R FP =R* P P = CpsiR's 5 + R s s p ) 


(7. 2. 2-2) 


These expressions provide the information necessary to define the fixed-frame con- 
straint. In addition, the copy air mass constraint can be defined at this time. In those 
cases where the superordinate frame is not the same as the parent frame, it is necessary to 
create copies of the independent and element air nodes in each of the subsystems leading 
to their nearest common ancestor. These additional air nodes are also connected together 
using copy air mass constraints. 

If the model containing the air mass connectivity constraint is to be used in a steady- 
state problem, the two cyclic air node degrees of freedom are meaningless. Therefore, they 
must be eliminated by defining two prescribed constraints in the superordinate subsystem. 
If, on the other hand, the model is to be used in an eigensolution, the gradient degree of 
freedom is meaningless. A single prescribed constraint is then defined in the superordinate 
subsystem. 
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7.2.3. Periodic Structure 


The purpose of the periodic structure constraint is to provide a simple means for 
creating an axially symmetric structure. This is accomplished by replicating a single branch 
of the model at equal azimuth angles about an axis of symmetry. For this constraint, the 
parent subsystem represents the assembled periodic structure and the child subsystem 
represents a single component. The periodic structure constraint consists of one or more 
of the following: a periodic node demotion constraint, a periodic generalized coordinate 
constraint, and a periodic air mass constraint. Note that the periodic frame constraint 
must be defined separately. 

When there are nodes in the component, periodic node demotion constraints are 
needed to transform them into the assembled structure. If the independent node corre- 
sponding to a dependent node (in the component) does not exist in the parent subsystem, 
a string of images of the independent node are created in the intervening subsystems and 
chained together with structural node demotion constraints. Similarly, if the dependent 
node does not exist in the child subsystem, a string of images of that node are created 
and chained together. Since the independent node (or its image) now exists in the parent 
subsystem and the dependent node (or its image) exists in the child subsystem, a periodic 
node demotion constraint can be defined. 

One or more periodic generalized coordinate constraints are needed if there are gen- 
eralized coordinates in the child subsystem. Similarly, one or more periodic air mass 
constraints are needed if there are air nodes in the component. A process identical to that 
used to connect structural nodes is used if the dependent and independent air nodes are 
not in the child and parent subsystems, respectively. 

7.2.4. Rigid-body Connection 

The purpose of the rigid-body connection constraint is to provide a simple means for 
connecting two nodes together rigidly. It is actually a special case of the screw constraint 
in which the translation and rotation degrees of freedom are both locked. 

7.2.5. Rigid-body Mass Connectivity 

The purpose of the rigid-body mass connectivity constraint is to provide a means 
for attaching a rigid-body mass element to a GRASP model. The element subsystem 
consists of a frame oi reference and a center-of-mass node, both of which must be connected 
to their counterparts in the existing portion of the model. Therefore, the rigid-body 
mass connectivity constraint is made up of a fixed-frame constraint and a structural node 
demotion constraint. 

In the definition of the rigid-body mass connectivity constraint, the position arid 
orientation of the dependent, element center-of-mass node C relative to an independent, 
existing node I is provided. The position and orientation of the independent node relative 
to its subsystem frame of reference (the superordinate frame 5) is also known. Then, the 
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position and orientation of the center-of-mass node relative to the superordinate frame can 
be written. 


tCS _ gClQlS 


R c s s ^C SI Rf 1 + Rs s 


(7.2. 5-1) 


After locating the parent subsystem of the element subsystem in the system organization 
tree, the position and orientation of the parent frame relative to the superordinate, frame 
can be calculated. Since the element frame and the element center-of-mass node are 
coincident, the position and orientation of the element frame relative to its parent frame 
can then be determined. 


C FP =C cp = C cs C sp 
Rp P =R C P P = C ps {R c s s + R S S P ) 


(7. 2. 5-2) 


With this information, the fixed-frame constraint can be defined. In addition, all of the 
position and orientation information is available to define the structural node demotion 
constraint for the center-of-mass node. In those cases where the superordinate frame is 
not the same as the parent frame, it is necessary to create copies of the independent and 
element center-of-mass nodes in each of the subsystems leading to their nearest common 
ancestor. These nodes are chained together using a series of structural node demotion 
constraints. 


7.2.6. Rotating Structure 

The purpose of the rotating structure constraint is to provide a simple means for allow- 
ing one subsystem to rotate relative to another. For this constraint, the parent subsystem 
represents the nonrotating structure, while the child subsystem represents a rotating struc- 
ture. The rotating structure constraint consists of one or more of the following: a rotating 
node demotion constraint, a rotating generalized coordinate constraint, and a rotating air 
mass constraint. Note that the rotating frame constraint must be defined separately. 

When there are nodes in the rotating subsystem, rotating node demotion constraints 
are needed to transform them into the nonrotating subsystem. If the independent node 
corresponding to a dependent node (in the rotating subsystem) does not exist in the parent 
subsystem, a string of images of the independent node are created in the intervening 
subsystems and chained together with structural node demotion constraints. Similarly, if 
the dependent node does not exist in the rotating subsystem, a string of images of that, 
node are created and chained together. Since the independent node (or its image) now 
exists in the parent subsystem and the dependent node (or its image) exists in the child 
subsystem, a rotating node demotion constraint can be defined. 

One or more rotating generalized coordinate constraints are needed if there are gen- 
eralized coordinates in the child subsystem. Similarly, one or more rotating air mass 
constraints are needed if there are air nodes in the component. A process identical to that 
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used to connect structural nodes is used if the dependent and independent air nodes are 
not in the child and parent subsystems, respectively. 
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8. ELEMENTS 


The GRASP element library currently contains three elements, the aeroelastic beam, 
the air mass, and the rigid-body mass. 


8.1. Rigid-Body Mass 

In GRASP, rigid bodies are modeled as being influenced only by inertial and gravita- 
tional forces. 

For the purposes of modeling the motion of a rigid body in an inertial and (possiblj ) 
gravitational field, consider a rigid-body mass element B that has an inertia dyadic /. 

Steady- State. The rigid-body mass element (fig. 14) has a body-fixed node JV and a 
frame of reference P. Node AT is initially coincident with the deformed frame F (P — 

and C F ' N — A). The virtual work at the deformed node N‘ is 

tfW = F N ' ■ 1 SR NI + M n ' • H n>i (8.1-1) 


from which nodal forces and moments can be derived. The nodal virtual displacment and 
rotation variables for this element are 8R ^ y ^ and St l>^ , respectively. 

The inertial virtual displacement and rotation of the deformed node N' are 


'iB? 1 = W ' 1 + s± Fl1 X r n ' n + n sr n ’ n 

S^' 1 =Stl> N ' N + 6± nf ' +H F>I 


( 8 . 1 - 2 ) 


The force acting on the body at N' is 


F N> — —mA N 1 + mg 


where the inertial acceleration of N' is 

a n i = a f '' + n 1 " 1 x (n F>I X r n ' n ) 


(8.1-3) 


(8.1-4) 


Substituting equation (8.1—4) into equation (8.1—3), and transforming from the body-fixed 
(N) coordinate system into the deformed-frame (F 1 ) system, 


jpN' / aF' 1 C)F'If)F'I jjN'N\ 

=m{g /v ~ A n - U N U N ll N ) 
( . F'l nF'r^F'I nN'N\ 

=77i\^gpi — Apt — s Ipt Upt 1 i/v ) 


(8.1-5) 
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Figure 14. Rigid-body mass element 


The moment acting on the body at N* is 

M n ‘ = H nI1 xQ n>i 

where the inertial angular momentum at N' is 

h "' 1 = i • n N>I 

and the inertial angular velocity at N 1 is 

qN'I =a N'N + n NF> +Q F'I 



The frame force and moment components can be derived in a similar manner. If the 
frame virtual displacements and rotations are 8R F , F and 8ipp, F , respectively, 


pF' 
r F , 

= C F ' N 

Flf 


II 


M F ', 

II 

+ R%' N Fif 


II 

►*** 

** 

&N' nN’N 
— r N U N 


( 8 . 1 - 10 ) 


Dynamics. Assuming that the rigid-body mass node is perturbed from its steady- 
state position (fig. 14), the virtual work at the node may be expressed as 


m = i>R%" lT Fif + s^S" iT mS" 


( 8 . 1 - 11 ) 


where the force and moment are 


the angular momentum is 


Fjf = - mAn' 1 + mg N 


Mf = - iNaf 1 + nf 1 ^" 1 


Hit" 1 = i N n»" r 


and the inertial angular velocity is 


(8.1-12) 


(8.1-13) 


jf'^nfr+nS'^' + c^np 


(8.1-14) 


The virtual displacement is then 

6RZ 1 = SRf"> F ' +*Rn 


n"i _ x a p n-n' + + r n-n + r n'n ] 


Lyy — VAt F 
and the virtual rotation is 


(8.1-15) 


sts" 1 = +w%" n ' 


(8.1-16) 
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The inertial acceleration of the node is 


./v”/ 

l /v 


. F"F' 


■F"F' 


=a f f , 1 + n*F, I n F p, I R% N + r f „' -r%' n 0 f „ +r n 


N n N' 


\R%' N n F F : i + (n' F ,'R%") +n t F , I R%”}0 F „ + 

N N +n£ I n£ 1 R$” F '+ 

[i£ ; - (n£'<"fn£' - npAg" n F F : i }e F F :' F '+ 

f)F' I fjf 1 / nN" N' 

M pf i L pi 1 Lyy 


\F‘ 1 oN'/vr 


- F'l nN'N 1 


■F"F' 


2QP’, 1 R n 


and the components of the gravitational acceleration are 


" /IF n F* 
9N = 9F‘ + 9F'0f" 


■h 2J7 


• F" F' 
F'/jS 
P> Up" 


( 8 . 1 - 17 ) 

( 8 . 1 - 18 ) 


The components of the inertia dyadic in the nodal basis can be expressed in matrix 
form as 


~N N 


i N = { a + 0 n )c NN 'i N ..c n ”{A~e N ) 


N'N/ 


-N"N' 


( 8 . 1 - 19 ) 


Finally, the angular acceleration is 


♦♦ F n F l •'N n N t • F n F f • N n N f 

ofyv ^ — 0^11 + o N 4" o j? lf -(- ri^^Of 


& pi V pti 


L F* u N 


( 8 . 1 - 20 ) 


The force and moment can then be obtained from the substitution of equations (8.1- 
12) through (8.1-20) into (8.1-11). When the virtual work is calculated, it consists of the 
same steady-state residuals {Q} as were obtained in the previous section, in addition to 
the virtual work associated with the coefficient matrices [M], [C], and [ K ]. 





f 

r •« F n F* ' 


' . p n p f > 





6r f : f ' ) 

T 


R pn 

.. p n p 1 


Rpn 
• F" F* 


( pF"F' >| 
Jlpn 


= < 

syy 

6R»" N ' 

UiT"' J 

► < 

[MJ. 

e F n 

» N"N‘ 

r n 

■iN" N' 

► + [C] < 

e F „ 

•„ N"N' 
Rw 

■_N"N' 

■ + [K] < 

dF" F f 
v pn 

niV"iV' 

k n 

nN"N' 

K V N > 

► - { Q } 




. > 


. ® N 



j 


( 8 . 1 - 21 ) 
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where the M coefficient matrix is defined to be 


SR F f : F ' row: 

.. f"F’ 

Rp,, column: 
•■F" F' 

0 pti column: 

•i N" N' 

R n column: 


•■N"N 

0 N column: 

c ; F" F' 
d\pp„ row: 

•• F"F' 

Rp,, column: 


”F" F' 

0 F n column: 


• t N"N' 

R n column: 


~N"N 

0 N column: 

6R%" N ' row: 

F U F' 

R F , t column: 
»F"F' 

0 F ,, column: 

.. l\?” N' 

R n column 


-TV" TV 

0 N column: 


rn 

- mR%' N 

m 

0 

mR%' N 

- mRfl' N Rf!i' N 4 - C NN> In"C n ' n 
mR%‘ N 

C NN> In"C n ' n 

m 

nN'N 

— ITlxl^ 
m 

0 


(8.1-22a) 


(8.1-22b) 


(8.1-22c) 
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row: 


•; F"F' 

R F ,, column: 0 


9 F n column: C NN Iw<C N N 


.. /V" N* 

R n column: 0 


(8.1 22d) 


el N column: C NN 'l N »C N ' N 

the C coefficient matrix is defined to be 


c r>F" F’ 
cl tp,, row: 


F" F* 


R f ,, column: 2mf!^ 


FI 


. F” F* 

Oirn column: 


• yv"yv' 


- 2mfl£, I R%' N 


R n column: 2mfl£, i 


(8.1-23a) 


iN n N 

9 n column: 0 


c t F" F 
ojp F u row: 


. F " F > 


R F , f column: 2mi?^ IW U F , 1 


•F"F' 


9 F n column: - 2mR%' N n£ T R$ N - (C NN J N .C iV ' N fi^ / )" + 

ryNN'r /liV'iVn F'l , fsF f t x-r/V'/V 

L> M F i -j- 11 F t C 1 /V” C' 


(8.1-23b) 


N" N’ 


R n column: 2mR% N Qp, 1 


■ N" N 


0 N column: Q^ 1 C NN ' I N ,.C"' N - (C NN ' I N ..C N N n£' )' + 


N'N tr>NN'- 


iN'NnF' I\ 


c nn 'i n „c"' " nz, 


N'NnF'l 
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row : 


6R 


N"N' 

N 


F" F' 

Rp,, column: 
lF" F' 

B F „ column: 

•_ N" N' 

R n column: 


2mQp, 1 

-2 m(l£’RK' N 
2mQ F < 1 


(8.1-23c) 


•_N"N 

9 n column: 0 


6^^" N row: 


. f u F r 

R pn column: 0 


F . x-iAT/V' r s-iN'N r\F f I , rN/V' T piV'N 

0^,, column: C ^ I^nC ttpi + ^ 


(8.1 23d) 


N"N‘ 

R n column: 0 


‘.N"N A.F'IsiNN't /~iN’N /^iNN't riJV'Nn F' l , 

0 N column: Sljp,C nri In»C‘ - {C In "C ) + 


,NN' 


T piN' N f*F' I 
i/V"^ J Lpt 


and the K coefficient matrix is defined to be 
6Rp„ F row: 


bF"F' 

xlp/f 

column: 

mU'p'Jn'p', 1 


nF n F* 

vpn 

column: 

mlAp' 1 — gF' — ((lp> 1 R% N ) ^F' l — 

n$, I R%' N n$. 1 ] 

(8.1-24a) 


column: 



qN u N 

V N 

column: 

0 
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row: 


6 - 0 


F" 


R F „ F column: - m[Ap, J - gF< + (Clp/ilp, 1 R^ N ) - 


i FI 


-F'lnF'l r>N'N\- 






0 F „ F column: mR $ fj£! 7 - 

nf ^jy^ng 7 ] + nS 7 c JVJV 7j V »c iV,N nS J - 

(C NN 'I N .,C N ‘ N nj; 7 )'n£ 7 


column: - m[ip 7 - j F . + (n£ 7 np 7 J*jf ")‘ 


j$'"tt£ 7 ii£ 7 ] 


0%" N ' column: n^C^^C^n# 7 - 

fyF' I( r NN ' t riN'NaF' I\' 


c DATIN' 

o/iyy row: 


Rfit F column: milp^Qp, 


F l I f\F* I 


0p.. F column: m[Ap, J - gp> - {^lp, r R%' N ) f}£i 7 - 

QF'/nN'NoF'/i 
it pt Jttjsjf S t pt j 


IF' I 


R% N column: mClp, 1 Up, 1 


pt ilp t 


6*% N ' row: 


N column: 0 


R f , f column: 0 


0p,, F column: flp^C^ 1 ^ I^„C N ^^tlp, 1 — 

(f^NN'r nN'NnF'/rnF'/ 
V*-' iTV"^ ilpi ) lip, 

Rtf N column: 0 

9% ^column: Cl F , r C NN Ijv"C n N {lp , T — 

p,F' I (r ,N'N T pyN'NpfF' I\- 

ilpi IpfUKs itpt ) 


(8.1-24b) 


(8.1-24c) 


(8.1-24d) 


82 



8.2. Air Mass 


The air mass element models the momentum flow of air through a helicopter rotor 
disk. For this element, the rotor is assumed to be an actuator disk, and the flow field a 
cylindrical region surrounding the disk (fig. 15). The state vector for the air mass element 
is made up of the generalized coordinates for a single air node. 



Steady-State. Consider the air flowing steadily through a rotor. Reference 22 shows 
that the thrust dT acting on a differential annulus of the rotor (fig. 16) is related to the 
induced velocity v via a momentum balance such that 

dT = 4np a rv \V + v| dr (8-2— 1) 


where r is the radial coordinate of the rotor and V is the velocity of the rotor relative to 
still air (V is positive when the rotor is moving in the positive direction). The use of 
the absolute value of the sum of the velocities V + v assures that the differential thrust dT 
has the proper sign under all operating conditions. Integrating, the total rotor thrust is 


f R 

T = 4np a / v \V + v\r dr 


( 8 . 2 - 2 ) 
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Figure 16. Air mass element differential annulus. 


The virtual work done by the thrust on the air is 

>R 


6W = 4np a j v8R \V + v\r dr 


(8.2-3) 


where 6R is the virtual displacement of the air. The right-hand side of the equation for 
the virtual work can be discretized by letting v = U ^ + 7 j^r and 8R — 6Pf + v6<f>f r . Then 


6W =47 rp 


•S> 


+ 7i*r r ) V 4- U A + 7 i r r (6Pf + r6(f>f r )r dr 


~6P A 4irp 


■ f i(fi 


A + Tit 7 ’) ^ + U* + l\r r r ^ r + 


A , --4. 


(8.2-4) 




/•* 

l /1 r 47r P« J (Ci 


• 4 +7 1 V)|V r + £^+7iV r 2 dr 


Ft A 


Note that while the coefficient of in equation (8.2-4) is equal to the rotor thrust, the 
coefficient of 8<j)^ r has the dimensions of moment but no clear physical significance. 

The contributions to the bV\? (applied loads) side of equation (8.2-4) are determined 
from blade element theory, and are obtained by summing the contributions from each of 
the aeroelastic beam elements that make up the rotor. 

Dynamic . Simple models for the induced inflow dynamics, such as the one introduced 
in reference 23, have been shown to improve the accuracy of mathematical models of heli- 
copter rotor dynamics. The velocity of the air mass is idealized as consisting of a spatially 
and temporally uniform freestream velocity V y which is augmented within a cylindrical 
region by the steady-state inflow Uf induced by the rotor steady-state thrust, and by the 
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infinitesimal dynamic perturbations to the inflow induced by dynamic perturbations to the 
thrust, roll moment, and pitch moment of the rotor. 


For a differential annulus 
the momentum balance can be 
equations. 

tfW - 


of a rotor disk through which air is flowing unsteadily, 
expressed as a system of first-order, integro-differential 



2 7T 

2p a v \V + v\8P difr dr+ 


p a v8PdV c ff 


(8.2-5) 


^eff 


In order to intermix the air mass terms with the structural generalized coordinates in 
a single set of second-order equations, the perturbed air mass generalized velocities are 
expressed as the time derivatives of generalized coordinates. 


• A ■ A \A 

V = v£ + 7j t r + P\ - 12 rsin V’ + <i>13 r COS ^ 


( 8 . 2 - 6 ) 


where P, is the vertical component of the perturbation of the induced inflow velocity 
1 . a . a 

component at the center of flow, <^ 12 and 0 13 are the flow gradients at the center of flow 
in the x 2 and x 3 directions, respectively, and i/> is the azimuthal coordinate of the rotor, 
measured as a right handed rotation about the xi axis from the x 3 axis. The flow 7 direction 
is assumed to be positive along the X[ axis. 

In addition, virtual displacement of the air inside the cylindrical flow field is assumed 
to be 

6P = 6P* - 6(/>f 2 rsinV’ + cosip (8.2-7) 

where 8P^ is the vertical virtual displacement of the air at the center of flow, and 
and 8<j>f 3 are the cyclic virtual displacement components at the center of flow. 

Now, consider the expression v|I r + u|, where v — v + u(^). In seeking the linearized 
perturbation of such an expression, if V + v = 0 then v\V + v| = (v + v)|v|. Since there is 
no linear contribution in this expression, it may be assumed to be zero. Now, define 

i +1 for a > 0 

0 for a — 0 (8.2-8) 

— 1 for a < 0 

Then, 

v |V + ti| = [(F + n)n -f ( V + 2v)-i;]sgn(F + v) (8.2-9) 
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Since only the linear perturbation dynamics are pertinent to this problem, the contri- 
bution of the change in momentum per unit area term from equation (8.2-5) is 


» R. i*2 7T 


II 2 p a r(V + 2v)sgn(F + v)(P 1 - c^ 12 rsin ip + «/» 13 rcos rp){ 6 Pf + 
64 >i 2 rsm + 6 <pf 3 r cos ip) dip dr — 

J 4tt p a r{V + 2v)sgn(F + v)[P 1 SP + —{4> l2 ^\2 + 4> 13 ^ 13 )] dr = 

•_ A f R 

P, SP J 4 47Tp 0 J (V -I- 2v)sgn(V + v)r dr+ 

2tt p a {<j> 12 6<P? 2 + iiMii) J ( V + 2v)sgn(V r + v)r 3 dr 


( 8 . 2 - 10 ) 


The contribution of the volume term from 
inertia effect as calculated in reference 23. 


equation (8.2-5) is the virtual mass-virtual 


8 Pa 

3 


{R 3 


€ 3 )PUPi A 


16pa 

45 


{R 5 


e5 )(^12^^12 


+ <P 


13^^13) 


(8.2-11) 


From this development, the coefficient matrix for the generalized accelerations may 
be defined to be 


M = 


8 Pa R 3 

3 


1 


0 

0 


0 0 
TE R2 ( 1-T&) 0 

0 A^(l-rfr)J 


( 8 . 2 - 12 ) 


and the coefficient matrix for the generalized velocities may be defined to be 


C = 2np a 


2 J f R gr dr 
0 
0 


0 

J e R gr 3 dr 
0 


0 

0 

/* gr 3 dr 


(8.2-13) 


where g — ( V + 2u)sgn(T -f v) and v — U j 4 + ry^. 

To eliminate all periodic coefficients in the equations of motion, and to assure the 
existence of a steady-state solution, the air mass element degrees of freedom must be 
inertial. In addition, the flow direction must be coincident with the steady-state spin axis 
of the rotor and the gravity vector, if gravity is included in the model. 


86 



8.3. Aeroelastic Beam 

The aeroelastic beam element is designed to model a beam undergoing small strains 
and large rotations, and for which shear deformation and warping rigidity may be ignored. 
A model of this type is developed in reference 24, which formulates the nonlinear beam 
kinematics and applies them to the dynamic analysis of a pretwisted, rotating beam el- 
ement. The kinematic relations that describe the orientation of the cross section during 
deformation are simplified by systematically ignoring the extensional strain compared to 
unity The only restriction on the magnitudes of the orientation angles used in describing 
the cross section orientation is that they remain less than 90° . All influences of warp other 
than warping rigidity are retained. The beam cross section is not allowed to deform in its 
own plane. The static equations from reference 24 are used without simplification; the dy- 
namical equations are linearized relative to static equilibrium. One noteworthy feature of 
the derivation of the equations in reference 24 is that the common practice of using an or- 
dering scheme has been abandoned. Thus, all liiglier-order terms (within the assumptions 

above) are retained. 

In the following sections, the details of the derivation of the equations for the aeroelas- 
tic beam element are presented. First, a synopsis of the basis under which the governing 
equations of the beam are derived is given. Next, the equations of motion for the beam 
element are derived in terms of the frame, air, bending, extension, and torsion degrees of 
freedom. These equations include contributions from beam elasticity, inertial and gravi- 
tational forces, and aerodynamic forces. Then, the discretization of the beam degrees of 
freedom is presented to show how the beam displacements are transformed into the beam 
generalized coordinates. The final two sections describe the transformation from root and 
tip node degrees of freedom to beam generalized coordinates, and the transformation from 
beam generalized forces to root and tip node forces and moments. 

8.3.1. Basis of the Governing Equations 

Consider the beam element shown in figure 17. The element frame is denoted by 
F, and the root and tip nodes are denoted by R and T, respectively. The addition of 
primes and double-primes signifies the static and perturbed dynamic states, respectively. 
It should be noted that F" and R are coincident with each other and that their coordinates 
line up with the principal axes of the root end of the undeformed beam element with the 

undeformed beam lying along 6? . Similarly, T is at the tip of the mdeformed beam 
element and its coordinate directions lie along the principal axes for the tip cross section. 
The air nude, denoted by A, must be included in the problem so that the influence of 
aerodynamic forces on the air node generalized forces can be determined and so that the 
influence of perturbations of the air node generalized coordinates can be determined for the 
generalized coordinates of both the beam and the air node. The position and orientation 
of A are inertially fixed. 
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Figure 17. Aeroelastic beam element (undeformed with pretwist). 

Interior displacements of the beam are represented by four functions of the axial 
coordinate x%: U{ and 0^. Bending is described by ui and 1 * 2 , axial displacement by u 3 , 
and torsion by These functions are discretized in terms of standard cubic and linear 
polynomials so that the generalized coordinates at the root and tip of the beam can be 
related to the nodal displacements and rotations. In addition, however, there are also 
generalized coordinates, called internal degrees of freedom, associated with higher-order 
polynomials. 


8.3.2. Beam Elasticity 

The derivation of the equations to calculate contributions of the elastic deformations 
of a straight, pretwisted beam follows the derivation presented in reference 24. 

Steady-State . The elastic beam equations for a beam in equilibrium are derived from 
the variation of the strain energy 


SU — f (Gej a Se^ a + £633^633 ) dx$ (8. 3. 2-1) 

J 0 

where 

^31 — (^1 — £2) (^3 “ 0 *) 

f32 = (-*2 + 6 ) («3 - 0 ) (8. 3. 2-2) 

f 33 =f33 + ^2«1 — £l«2 + ~ (£l 2 + ^2 2 ) ( k 3 — O') 2 + (^2^1 ~ £ 1 ^ 2 ) («3 “ O') O' 
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where 0 (^ 3 ) is the pretwist angle (fig. 17), with 0(0) — 0, and ( )' — d( )/dx 3 . The 
generalized strains are 


^33 =S' - 1 


s' 2 =ui 2 + «' 2 + (l+^) 2 

ki =(Cj2tt'i' - (8. 3. 2-3) 

K-2 —{C22 u 'i ~ C21 u 2)/C33 
( Cs\ 2 Cs2U 1 p 

«3 =^3 - ( - _ g a ~ + C 31«2 J / G 33 

where C = C P F ' , the direction cosines of local principal axes relative to the static frame 
orientation. The elements of C may be expressed in terms of Tait-Bryan orientation angles 
(orientation angles of type body-three: 1-2-3) as 

=C 2 C 3 

C 12 =s 2 ci + 51^2^3 (8.3.2-4a) 

C13 =S 3 5 l - C1S2C3 


C 2 ] — — C2$3 
C22 =C3Cl - ai3 2 S3 
C 23 =C 3 S] -|-CjS253 

(8.3.2-4b) 

G 31 — u \ 

C32 =u> 2 

C 33 =(1 — Uj — «2 ) 1 
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where 


5 j — Uj (1 ^] ) 




s 2 =«1 


5 3 = sin 0 3 


2^5 


C, =(1- V) 
C 2 =( 1 - s 2 2 )\ 


c 3 = cos 0 3 

After integrating over the cross-sectional area, the variation of the strain 
obtained in terms of the stress resultants F 3 , M j, M 2 , and M 3 . 


s: 


where 


6 U — I 6 s f *4* Mj Sk] 4- M2 M 3 6 k 3 ) dx 3 


F 3 = Eq€ 3 s 4~ — E\K 2 4“ ~ r 3 2 + DqQ* T$ 

& 

Rn Yl 2 

Mi =E 2 e 33 +hKi + -?-±- + D 2 0't 3 

B r 2 

M 2 = — E 3 ( 33 I 2 k 2 — 


, B 3 t 3 2 3 D 3 d' 

M 3 = I J + I 3 C 33 + B 2 k\ — B\k 2 H — 1 — — 


r 3 + D a 9' 1 r 3 + 


1 


(£> 0 e 3 3 + D 2 ki - D\k 2 )6’ 


(8. 3. 2-5) 


energy is 
(8. 3.2-6) 


(8. 3. 2-7) 
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where r 3 = k 3 - O' and the section integrals are defined as 

E 0 = J J EdA 
E[ — J j E(idA 

E 2 = JJ E(, 2 dA 
Ji = JJ E(, 2 2 dA 
1 2 = 


//*>*" 


I 3 + ^2 


J = JJ G[(Ai - 6) 2 + (>2 + ti) 2 ]dA 

Bi = fj £ 6 ( 6 2 + 6 2 )<^ 

B 2 = J J EMti 3 A ^2 2 )dA 
By = J J E{it 2 + h 2 ) 2 dA 

Do=JJ E{( 2 X 2 - ^X 2 )dA 

Di=JJ E^Xi - hX 2 )dA 

D 2 =JJ Ehib* 1 ~(iX 2 )dA 
D 3 = J J £( 6 2 + 6 2 )( 6 ^i - (iX 2 )dA 
D4 = J J E(( 2 Xi — £iX 2 ) 2 dA 


(8.3.2-7b) 


91 



Here, Ey is the axial rigidity; E j and E 2 are the first flexual moments about the local 
£j and £2 axes, respectively; Ij and I 2 are second moments ( bending rigidities ) about the 
local £1 and £2 axes, respectively; and J is the Saint- Venant torsional rigidity. 


The variations of generalized strains can be expressed in terms of the fundamental 
variables as 

, ds' 

6s' = - — Su t 


dui 


Q pc ' 9 k, ‘ 

fifci ~ bu^ + ^ e 3 ia K ab6$ 


(8. 3 . 2 - 8 ) 


dul 


and the variation of strain energy as 


bU 


/[( 


ds' dn \ ds' 

F 3 - ( 4- MilTT ) ^ u 'a + E 3 Q^T^ U 3 d" e a03Ma K 0^3 + 


du[ 


du ' 


M iir ±6u”+M 3 S 0 ' 3 


d K 


(8. 3 . 2 - 9 ) 


dx 3 


where 

ds' S 3 i + u' 
chi' s' 

dn 3 E \ 2 

d^ = C^ 3 
dn j C n 

d^ = ~c7 3 

dn 2 C22 

d < 

dn 2 _ C 21 

~~ cT 3 

9 k 3 — C 3 \ 2 C 32 

du" t-^33 ( 1 — C31 2 ) 

dtt 3 _ C 3 \ 

du C33 

9 k 1 _ u" / (^ 22^32 C12C3I \ U3 t^i)C3iC| 2 

9 u'j C33 2 \ 1 — C 3 J 2 C33 J C33 3 (l — C 31 2 ) 


( 8 . 3 . 2 - 10 a) 
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d*2 U| / C\ 2^32 C2 2C31 \ 1>2 C2] C31C32 

"c^ 2 \1-C 31 2 C7 33 / _ C 33 3 (l - C 3J 2 ) 


dn> _ -c„c, 2 u'; 7 + Cji 2 _ CjiJ) 

C 33 3 (1 -C 3 1 2 ) 


tt 2 , ( 1 ~ C32 2 ) 

c 33 3 


(8.3.2-10b) 


Oki 

du' 2 



(u"C 3 i + t ^2 ^32 ) 


Dynamic. Since the explicit, analytical derivation of the elastic stiffness matrix 
would be exceptionally tedious and lengthy, GRASP generates it numerically. This is 
accomplished by taking the Jacobian of the function that calculates the steady-state elastic 
loads. Because of the necessity of calculating an accurate stiffness matrix, the algorithm 
used to calculate the Jacobian uses a two-point central difference scheme plus a generalized 
formulation of Richardson extrapolation. 

8.3.3. Beam Inertial and Gravitational Forces 

The generalized forces resulting from motion of the aeroelastic beam relative to an 
inertial frame are also determined following reference 24. Warping dynamics are again 
ignored. The derivation is based on the work done by inertial and gravitational forces 
moving through a virtual displacement. The work is calculated by taking the scalar product 
of the gravity minus the acceleration of a generic point P (fig. 17)in the beam interior 
(ref. 24, eq. 32), with the virtual displacement of the same point (ref. 24, eq. 34), then 
integrating the result over the beam length. 

Steady- State. For a beam element in equilibrium, the virtual work is 
8W = / ( Su r Fp + Sifrp.Mp, )dx 3 + 

Jo (8. 3. 3-1) 

f t »t v 

Sup , I Fp dx 3 + Sxj ')]?< I (Mpt + Rp> F Fp, )dx 3 
Jo Jo 
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where 


Ff, —m(g F i A^, 1 ) - il$ I n$ I R% F ' 

M pi j — Tn 2 ^P '3 — 

AIpt2 —— ^1 Gp'3 -{- 

Mp, 3 —vtixG pt 2 — m2t?p'i 4" (*i — i2)^p # i^p f 2 
Rpif =m(x 3 6si + Ui ) + mjCji + m 2 C 2 i 


and 


Gf> =g F , - A*', 1 - QpnpRZr 

u =R% P 

6u =6R^, P 

64 > P ' t = 6 M + Pli s < = WpT 

Supi =6Rp, F 


The section integrals are 


SiftF 1 


- 6 ift F , F 


m = 



p,dA 


777.1 = 


// 


PtCidA 


( 8 . 3 . 3 - 2 ) 


( 8 . 3 . 3 - 3 ) 


77l 2 



*1 



*2 



( 8 . 3 . 3 - 4 ) 
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Dynamic. For a linearized perturbation about the equilibrium solution, it is possible 
to express the equations of motion in a matrix format such that the virtual work per unit 
beam length is given by 


( bup"i y 

T 

r 

•• > 
UF"i 


t i ' 

UF"i 


f Up**i y 

> 




0 F" i 


0F"i 


0F"i 


Sui(x 3 ) 

> < 

[M}< 

Ui(x 3 ) 

> +[C}< 

Ui{x 3) 

> + < 

Ui{x 3) 

► - {Q} 

tu' a (x 3 ) 



Cm 


Cm 


C(* 3) 


vS0 3 (x 3 ) ; 


! 

,^>M ; 




J 3 {X 3 ) - 

J 

/ r% 


where the components of the generalized force vector Q are the same as the static general- 
ized forces (see the previous section) and the coefficient matrices M,C, and K are defined 
on the following pages. The M coefficient matrix is defined as 


bup"i row: 

upn j column: mSij 

Opnj column: tijk(mRp F ) 

Uj column: m6 t j 

.. ^ ^pi prt dni 

u '{3 column: 

" p t p > 

0 3 column: 


(8.3.3-6a) 
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bipF"i row: 


0F"i column: 6ij(mRp R P 'k + 2 miCfe* ) 


iP' F' nP'F' 


m pP'P' dP'f' m irih 
Tnitpi t Kp t j — 


P' F* 


r p f :f + c, 


P’F‘ 


pP’F' 
K F f i 


Uj column: 


c tjk (mR 


P*F* 

F'k 


-j - 771 1 C 


P’F 

Ik 


(8.3.3-6b) 


2, column: C£T - C^f C* r ' )R<*. £ 


+ 


, r P’F‘ i ■ r>P'F‘ , ■ r P'F' & K 3 

lCu dO» + * Cai ^ + 3 ° 3 ‘ 


du 


0 


du^ 


0$ column: F C k{ F 


iP'FVP'F'n nP'F 1 


31 


ki 


)R$> 


+ 13C 


P'F* 
3 i 


bui row: 


up’fj column: 


mb ; 




0F"j column: 




Uj column: mb t j 


U{3 


column: 


p>p< 9 k\ 

tklrnrn m C ki — ■ 

du 0 


(8.3.3-6c) 


8 3 column: 


d' p* 

« k3lTTliC ki 
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du' row: 


.. t /^T P* F' 

u F"j column: eklm^m^kj 


du” 


k 1 /riP'F'^P'F' r P'F' r P'F‘ \pP'F’ dK k . 

8 F «j column: mi{C kj C lm - C,_, c fcm )it F , m + 

• & K l ^P'F' , & K 2 r,P'F' , _• 5k 3 ^p’F' 

+ ,J S" C » +,3 a< c >' 


” p'F ‘ 9 ki 

Uj column: efei m m m t/ fcj 


Ok i 5/t] . 9^2 <?/C2 • ^ K 3 <^3 

^column: i, ^ + «3 ^ 


$3 column: 


1*3 


d*3 


6 3 row: 


ujr»y column: 


eksimtC 


P'F’ 

kj 


"a 1 {n P'F' r P'F' r P'F' r P'F' sjjP'F' , 

column: ^.*(£3^ C kl - C 3 , C fcj + * 3 '-'jj 


tP'F' 


ttj column: ekil m iCkj 


, . #*3 

u^ column: 13 -jr-jj 

ottp 


63 column: 13 


(8.3.3-6d) 


(8.3.3~6e) 
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the C coefficient matrix is defined as 


Supi’i row: 


upfj column: 




0 F" j column: 26ijfl F , k (iiiRp , -f rn l^'[k F )' 

2Qp, I j (mRp f f + m k C ki F ) 


u 


pf r 

j column: — 2€ijkmflpi k 


t w?t v 9 k t 


~ 0 F'// r TP'FV p '*’' piP'F'riP'F* 

Up column, 2m*S£ F ,j (C mi G u — C ifei G m j ) ^77 


(8.3.3-7a) 


column: 2 m k Sl£f (C F - F ‘ C% F ' - C£ * C FF ' ) 

Uf>j column: - 26i } Q F i l k (mR F , k + mjCj^ f )+ 

2ft£f(mHgf + m*< F ') 

- rl r n» rl n< pf 

ffpnj column: — 2e l jjtmf! F ,j ii F ,j — 

9 -r ^ nF'//rP'F' dP'F' , r P'F' nP'r' \ 

**ijk m m**F'l ^F'* + C m* ^FW /” 

2 Cijfc *3 fc 


6lj)F"i row: 


Uj column: — 2£^Q£ f £(mi?£,£ -f m/C^ * )+ 


p' F* 


,P'F # 


2fl£,7(mfl£,'f + m*C , £ i ') 


P'F* 


Up column: 2e* Jm Tn m (ft£,fc£ l F - fl F . l n C ki F )Rp, F 


0 


2 n^(e tlJll C + eu,t 2 Cfr )C£ F '|^ 

du'p 


03 column: 2ej t omi(fijrifC'£ n F — * )Rpi^ n -f- 


F'/ riP’F' \nP' F' 


oo f'//j r’P'F'r’P'F' i r'P' f' r<P' b " \ 
ZM F ‘k\ l l l -'2k °lt t/ 2i ) 


iP 1 F 1 r'P' F' 


P'F‘ r'P' F‘ 


(8.3.3-7b) 
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6ui row: 


* o F' I 

UF"i column: — 2cij k mU F i k 


dp«j column: 26ijU F t k ( / tnR F , k + rajCj* F ) — 

2 n£<(mJ#r +m l c£ f '') 


F'/ 


Uj column: — 2 


(8.3.3-7c) 


l> , O 0 F'//^r'F , r -P'F' r P'F' r P'F'\ dK "' 

Up column: Zm k U F ,i{^ mi u fci - u mi ’ du'p 


0 3 column: 2m k n F ,[(Cf i F ^ kl * ~ ^3 1 * t ) 


rP'F’nP'F’ 


6u' a row: 


Ufnj column: 

0 0 F'// r ,P'F>P'F' r P’F' r P'F'\^m 

-2m k \l F n\ L 'mi '-'kl kj ^ ml ) q u „ 

0F"j column: 

- 2t klm rn m (n F F : i j cC F ' - <> T n C F ; F ' )R P F F ' l 


2si + «»i^r: F ')< F ' f| 

Uj column: 

o nP’t(nP'P'r p ' F ' r P'F , r P'F\ dK m 

- 2m fc SI F ,, (G mi C kl u mi ; 5u / f 

Up column: 

r-,# r r>f rpt . p* J7 f ^ 

- 2n Flm (ejfcf 2 riC 2m +Cfcji*2C lm )^ u „ 

6 3 column: 

r p'F '^ K i i^r pl f '^ K2 -) 

2n F , k {iiC 2k 


(8.3.3-7d) 
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60 3 row: 


upn j column: — [C F - F C F t F — C F t F C F j F ) 


d F "i column: - 2e fc3 /mi(n£^C’£] l F ' - ttp,l n C£- F ‘ )Rp,^ - 


90 F>I (i r ,p ' F ' r p ' F ' _ ; r ,b ” t 'r r t \ 


P'F'r'P'F' 


is column: - 2m k n F '.i(C 3 p ; F ‘ F ‘ - C F ' F '< F ') 


(8.3.3-7e) 


^ column: - 20^(1,^ ^ 


i F' l /■ /iP'F' 

du 


0 


0 


03 column: 0 

and the K coefficient matrix is defined as 
bu F " row: 

up" column: mClp.^p, 1 


0 F " column: rn[.A F , 1 — gp< -f R F , F ) 


n F : / n F : / ^ F : F 'i+ 


) - flp, I Upi I mF- 


- F' InF'l. 


(8.3.3-8a) 


u column: mflp^Up, 1 


0p> column: — Qpt^flp^rhpiC 


yF'P' 
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row: 


u f .i column: - m^A^, 1 - gF 1 + (ftp^&F' 1 Rp< F ) — R P > F (l F ,H ] 


(J7 p/Up^mF 1 ) + m/rfftp/Qp, 


- F'lnF'l 


9f« column: m{Rp, F {Ap, 1 — <7F')+ 

R£ F 'W£. I n£. I R£ F, )~ - n£, I n£, I R£ F '}}+ 

- IaF'I - i (nF' In F ' I nP' F' 
mpr\Api — gF' +{il'F' “F' n F' ) — 

nfI / nf! / A?! F ']+ 

Rp, F ^(Up, 1 Clp, 1 mp>) — fl F , I ilpt I rhF'} + 

riF'P'tTT o F'l c\F' I : ciF'/ \/~iP r F 1 

L/ 1 It pt J l pi 


n^ip’np , 1 ^ 1 


(8.3.3-8b) 


, r IF'/ - , / r\F’ I r\F' I r>P' F' nKf nf fnr 

u column: — m\Ap, — gp' + “f' R-f' ) Rp' U F' 1£ F' 


-P’F'nF ' InF'I] 


(Up.^p^mF 1 ) +ThF>Up, 1 Qp, 


~ F' InF’1 


6p> column: — R PF ilp^tlp^mpiC 

(m F .n£ , n£'R£ r ')'c r ' F '+ 

I T ~ TP* T — . D 1 T?* v P' 


T F'P' 


column: 


m iw (nf: / nF, / iif; F )‘C J 

[mF.(Apl f - gF ' )] C FP ' + rhp>(Ap, 1 - gF')C F P + 

c F ' p ' (Hp'Si'p, 1 - n% I i P >n F ', I ) T 


- pi | ■■ p 1 1 

Upn column: 


0 F ,< column: m[ip/ - i?F' 4- {Up, r (l F , T Rp. F ) — ^f'^F' 1 Rp' F ] + 

(8.3.3-8c) 


) — fljp # 1 fl p/thp 


~F'/nF'/ ~ 


u column: 


Opt column: — fljr/f l^JrhptC 


iF'P' 
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SipP* row: 


pf jp9 ^ ^ ril r ~ rtf r 

up" column: C mp<!]p Hp# 


Qpn column: C p F 77ipt(>lp # 1 — gp > ) + C l * rhpi^flp, 1 {lp f J Rpi h ) 


\F' I 


P’F ' - 


F’lnF’I nP'F 


I C1# - 


+ (tf P .n£;' - n^;'tp>n|>7)c 


r r' r “ r f 


,F'/\/-.P'F' 


(8.3.3-8d) 


u column: C F F mpiilpi^p, 1 


where 


0 P < column: C p F rhp^Ap, 1 — gp> + (Qp, 1 Sip, 1 Rp,* ) ] C* F -f 


i F'i 


-F’lnF’l nP'FVl^F'P' 


~ F J r ~ Z7 1 * 


i/p»llp#^ — Ilp/tp'fip# 


F /; 


>F / 


Hp> —ipittpt? 


I pi = 


M 

0 

0 


^2 

0 


0 

0 

*3 


(8. 3. 3-9) 


m p/ 


mj 

m 2 

0 


In the foregoing matrices, m is the running mass per unit length, and m a is the first 
mass moment about the £ a axis. The last block row associated with Sx^p, p is used to 
obtain the terms associated with 8u* a and 86 3 by substitution from the equations 


= hiMi + (^- + 


8k, 


8 2 K t 


- f - 

a< T au»au^ U/ * 


# 2 /c* 


Oi)Su a 


(8.3.3-10) 


and 


r 


p'» 


P' _ c -» 

- *«*» + 


(8.3.3-11) 


The 6iJ)p> row block matrices must then be pre- multiplied by R 1 and the Opi column block 
matrices must be post-multiplied by R, where 


R< 


8k{ 

du H 

a 


-^t3 — ^t3 


(8.3.3-12) 
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The geometric stiffness matrix K° is added to the 6tj>p> rows and 9 P > columns where 


rrG 


n' d 2 Ki 

Mp,i du^ p 



AfP' &*** 

p,i du^d0 3 


(8.3.3-13) 


K?i=* 

The matrix K G comes from the last two terms in equation (8.3.3-10), which are commonly 
called the geometric stiffness terms. 

8.3.4. Aerodynamic Forces 

The aerodynamic forces acting on the aeroelastic beam element are determined from 
a quasi-steady adaptation of Greenberg’s thin-airfoil theory (ref. 29). Before the theory is 
discussed in detail, two new sets of axes must be introduced for the purposes of defining 
the directions in which the lift and drag forces and the pitching moment act. In figure 18, 

the Z axes are associated with the zero-lift line for the airfoil section with the vector b 2 

z 

along the zero-lift line toward the trailing edge. The vector b 3 is along the beam axis but 
in a direction such that a dextral rotation of the airfoil section about ^iis vector results 
in an increase in the angle of attack. Then, being a dextral system, ki turns out to be 
normal to the zero-lift line (and nominally in the direction of positive lift for the section). 

The other set of axes is the so-called wind axes W. For these axes the base vector b 3 
is identical to b 3 . The base vector Sf is located along the relative wind vector (in the 

- w 

direction of drag) and b j is in the direction of lift. 


Zero-lift Itne 
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The Z basis and the P (principal axes) basis convect with the blade cross section, and 


are related by the direction cosine array C zp = C z p = C 


Z P 


C ZP = 


f cos 9 Z 
-a sin 9 Z 
0 


sin 9 Z 0 

a cos 9 Z 0 
0 cr 


(8. 3.4-1) 


where a — +1 if a dextral rotation about b 3 results in an increase in the angle of attack 
and <t — — 1 if a dextral rotation about b 3 results in a decrease in the angle of attack. 
The wind basis W is related to the Z basis by 


C wz 


cos a — sin a 0 
sin a cos a 0 

0 0 1 


(8. 3. 4-2) 


where a is the angle of attack. Then, C wp = C w z C ZP . 

Point Q is the quarter-chord point of the cross section, about which the aerodynamic 

forces and pitching moment are calculated. The offset position of Q relative to the origin 

* z 

of the local principal axes P is • 

Consider the wind velocity vector at the perturbed position of the aerodynamic center 
Q". W Q is calculated by subtracting the inertial structural velocity at Q n (V® *) from 
the inertial air velocity at Q n (U® 7 ), where 


• A 


~, A A 


and 


UQ"r = _ {Ui + t ^a + p j + jfQA^ + 

v*" 1 =(q f>i + n F " F ' + n PF )R QP + (f 1 + q f " f ')r pf " + 
f "r pf " + u f,i r f " f ' + f 'r f " f ' + 'if 1 


(8. 3.4-3) 


The relative wind velocity components in the Z" basis are then 

+ n A << + <<)- 

l(c z " F "n F :'y 4 - (c z " F "d F ..y 4 - {c z " p "'C p ')}r¥- 

R p ;r - 


*Z" F* 


( C Z P 9frt >) C z * Uf n — C 


ii — [{C zn F n ) C Z n u F ,i — 


Z tl F n - 


iZ"F” - r'Z"F"\zF'I 

U F u — ^ V pi 


(8. 3. 4-4) 
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and the local air flow velocity gradient is 


G% 


Z" 12 


Bwf, 

9R^ 


cf/( 


zai R*!? + R’ii 

V 

c*fn zF ; + c z r^ + c z tf 


W “ Aa 7 A + C az ^ 2 + C az ^)+ 


(8. 3. 4-5) 


where 


= \fn¥ 2 + 


If the time derivatives in equations (8. 3.4-4) and (8.3.4 5) are replaced with variations, 
the relative virtual displacements and rotations of an element of air with respect to the 
structure are obtained. 

6S%' ti = - C Z " A (SP A + RVm, + R-V^n b)- 
l(C Z " F " SlpF" )~ + {C z " p "Si, P P'' P ')]i2R q Z 2- 
(c z ' ,F ' , 6^); j c^ F,, Rp k F " - c?; F "f>u,- 

C z " F "6u F .. } - (C z " f "6x!> f .,)~C z " f "R f f : f ' 

= - c z : a c z ; a h a 2 - c?: a c% a h a \ + 

Z (8. 3. 4-6) 

The relative wind velocity magnitude and components are time-dependent quantities. 
For the magnitude note that 


w 2 = {wf,: x Y + (wg: 2 r 


(8. 3. 4-7) 


for which the static part is 
and the dynamic part is 


T / 2 — nlrQ' i 2 


n 


Q \2 


0*7., r + W, 2 ) 


w ~ if 


(8. 3. 4-8) 


(8. 3.4-9) 


Likewise, the angle of attack is a time-dependent quantity. In the equations written 
below, it is necessary only to develop the static part and the linearized dynamic perturba- 
tion part. These quantities are easily determined from the definition of a. 


tan q = 


C, 


(8.3.4-10) 
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The static part is simply 


tan a — 


*$! 


while the dynamic part is 


a = 


ryQ' wQ" _ w® W®" 

vv Z* 2 vv Z" 1 vv Z'l Z" 1 


W 2 


The applied force is assumed to be 


W n W h 


F = C c b x + Vb 2 + C nc bf 


and the applied moment is 


M - Mbf 


The equations governing the aerodynamic force components are 
£c =^p a W 2 cci + ^p a c 2 WG ^,,, l2 


V=-p a W 2 cc d 


m Jpjv’c'cn - * 6 p„c 3 (wg %\ 2 + + j6%: n ) 


C —~p c 2 (W + -G^ 1 

*" nc ~ 4^ at \* v Z " i ^ ^ Z" 12/ 


(8.3.4-11) 


(8.3.4-12) 


(8.3.4-13) 

(8.3.4-14) 


(8.3.4-15) 


Now, all of the quantities that are needed to define the virtual work are available. 

6W= {-SS% , , i Fz"i+ST%' l3 M)dx 3 (8.3.4-16) 

Jo 

Steady- State. The static generalized forces can be removed from the expression for 
the virtual work and written in the form <5W = J 0 * 6q T Qdx 3, where 8q is 


8q = 


6P A 

< 6up"i 

6lpF"i 

6ui 

6u' a 

Ue 3 


\ 


> 


(8.3.4-17) 
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and the elements of Q are 


: 

C c C AW ' + vcf"" 


r (C'CtF + VC?y ) - MC “ A (RVc** + Rfrcg*) 

Hn ■ 

0 

Wii-- 

0 

SllF'i ■ 

C'Ctf** + PC 2 T f ' = F t A 

bxl>F'i '■ 

MC?;r +(R% F 'F A )i 

8ui : 

f a 

8u' a : 

c ?f(M - 

60 3 : 

C? 3 p {M - r q z p F a ) 


( 8 . 3 . 4 - 18 ) 
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Dynamic. After removing the steady-state contribution to the virtual work, the 
virtual work per unit of beam element length done by the aerodynamic forces and pitching 
moment can be put into the following form: 






This equation can then be rewritten in terms of aerodynamic M , C , and K matrices, where 


M =AFH 


C =AEH + AFG + BH (8.3.4-20) 

K =AEG + BG + D 
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The elements of A are 


hP A row: 


6</>f r row: 


row: 


Wn row: 


C c column: 

n \V' A 

- '-'n 

V column: 

riWA 
~ l -'21 

C nc column: 

n Z’A 

~ '-'11 

A4 column: 

0 

C c column: 

0 

V column: 

0 

C nc column: 

0 

M. column: 

0 

C c column: 

r,W'A rjQ' a 
— t 'll 11 A2 

V column: 

r \V A nQ' A 
~ ^21 IV A2 

C nc column: 

n Z’ A nQ'A 
“'-'ll ^.42 

A4 column: 

n Z'A r ,Z'A 

'-'11 '-'22 

C c column: 

r \V A nQ' A 

~ U 11 11 .43 

V column: 

n W A dQ’ A 
- *-'21 n A 3 


C nc column: - 

Ai column: Cf x A Cf 3 A 


(8.3.4-21a) 


(8.3.4-21b) 


(8.3.4-21c) 


(8.3.4-21(1) 
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Sufr"i row: 


C c column: — C™ F 

T> column: — F 

C nc column: - C F 

M. column: 0 

Simp'll j row: 

C c column: F R^, F 

V column: F /?£-, F 

C nc column: t,jkC?j F 


M. column: — C ^ F 

Sui row: 

C c column: — C ™ F 

T) column: — C™ F 

C nc column: — Cf x F 


SA. column: 0 


8u' 


row: 


C c column: 


r u' 

L 'n 


Z' fiZP dQP ® k 
'-'33 11 z 2 


d< 


V column: 

Cnc column: Cf 3 p R% P 

M column: — — — 3 C? 3 P 

du » 33 


(8.3.4-21e) 


(8.3.4-21f) 


(8.3.4-21g) 


(8.3.4-21h) 
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60 3 row; 


r 1 „ riW Z' r>ZP oQP 
L c column; Cjj o 33 n Z2 


t> column: 


(8.3.4-211) 


C nr column: C 33 p R% 2 


A4 column: — 


The elements of D are 


SPf row: 


. _ „ w ’ Z' * Z 1 

-W&cg*jc e ^ + »*&-) 


W Zl column: 


IV z 2 column: 


|1C«|= 

Wf,C,f' 1 (£ c S£fe-+P^fe-) 

\wg\* 


(8.3.4-22a) 


G z i 2 column: 0 


6(/>f r row: 


W z j column: 0 


U 7 22 column: 0 


(8.3.4-22b) 


G Z12 column: 0 


6<j> j 2 row: 


column: 


column: 


iwfl* 


w?, C *' * R** (O + 


|wf| ! 


(8.3.4-22c) 


G Z12 column: 0 
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6<f>f j row: 


W% } coliimn: 


|tv«p 


1 i7j column: 


|wfl 2 


G < z \2 column: 0 


Supni row: 


Wg j column 


+ Z>g%~) 

|W|| 2 


'O iVg, C* a (C c '»> + £> ) 

W$ 2 column: V ^ — - 


\W?\ 2 


G Q Z12 column: 0 


Si/jpi'i row: 


W Z1 column: 


e t C z ' F ' R^ F ' ( C. ac n | T) ^£m ~i 

e tjk>* Z2 L/ lj n Fii e \L'c Qa + v g Q ) 


|W?I ! 


W% column: 


\ W %\ 2 


G z 12 column: 0 


£u, row: 


W % | column: 


lVg i C ; f^(£,^fe-+P^_) 

IWfl* 




oc!? 


11’^ column: 


|W«P 


Gj 12 column: 0 


(8.3.4-22d) 


(8.3.4-22e) 


(8.3.4-22f) 


(8.3.4-22g) 
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6 u' a row: 


Wg j column: 



Wg 2 column: 


firQ c zp nQP i 

iwfp 


IX- 


0<* 


+ v 


ec 


Z3L- 

d& 


Gz j2 column: 0 


8 O 3 row: 


W% j column: 




I V' 2 r 


dot 


+ V 


ec 


2X- 


ddt 


iwfl 1 


Wg 2 column: 




tr'*' 


(?J 12 column: 0 


(8.3.4 22h) 


(8.3.4-22i) 
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The elements of D are 


SP A row: 

P A column: 0 

< p A r column: 0 

<j>f 2 column: 0 

<j> A 3 cohimn: 0 

ujrttj column: 0 

0 F "j column: t tlj (C c C™' F ' + VC^' F ')C^ A 

Uj column: 

Up column: 

#3 column: 




W' P‘ \r>F' 




'm(C'C%r + VC? k - W 


W'P'\r,P'A 


(8.3.4-23a) 
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6<pf r row: 


jPj 4 column: 0 

<j>f r column: 0 

column: 0 

</>^ 3 column: 0 

Ufu column: 0 

6p"j column: 0 

Uj column: 0 

u'p column: 0 

9 3 column: 0 


(8.3.4-23b) 
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6(f) f 2 row: 


Pf 4 column: 0 

<f>^ r column: 0 

<j )^2 column: 0 

0^3 column: 0 


UF”j column: - (£ c C^ ,>l + VC% A )C? 2 


iW A\/~tF' A 


t a - /V p ' 

it 


@F" j column: - e iW [(£ c C^' + PC 2 w m ' z ')(C^C^^5?; 

siZ' F' f-iF' A f,Q' * l\, 

°mJ L '*l if /l2 1 + 

a jrr’Z' f 1 r-F' A/^z' * r'Z' Ar'F'Ar'Z'F' \i 
■ M { L 'U U *1 '-'22 — '-'11 '-'J2 *-'2fc JJ 

Uj column: - (£ C C*{" >1 + 

*i column: e klm [(C e C^ z ' + PC 2 T' * )(<?,?' A c£ A C 2 ^ Hjf - 

KA i n ZP r P' A n Z' A r Z' AcP‘ A r ZP\\^ K ^ 
7V H C lk C /1 °22 ~ L/ U L, k2 C 2J 

6, column: JC^CjT + PC^HC^C^'C^gf- 

cSrcK A s^)+ 

jt j ( Z P /-i P 1 A Z‘ A si Z' A (~iP' AsiZP \\ 

• /Vl l U lfc U f 1 U 22 ~ °11 L 'k 2 C 2 1 1 J 


(8.3.4-23c) 
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6<j)f : 3 row: 


P* column: 0 

<j)y r column: 0 

<j>^2 column: 0 

0^3 column: 0 


U 


ij column: — (C c C\\ A + T>C^\ A )C P 3 A 


Opi'j column: - ejkt[(C c C^ m z + VC™ m z ){C^ A C^ A Rfh. 

cSfc[; A R% A ) + 

m(c z ; f 'c[; a c z ; a - cfsc^cfS)} 

uj column: — (£ c C” A + VC]\ A )C P 3 A 

u' column: e t/m [(£ c C t T Z ' + VC]]" Z ' ){C Z ' A C& A C ZP R Q Z P 


r ZP r P' A j)Q' As , 
c /i u az 1 + 


M{C zf C?; a C z ; a - C Z [ A C P 3 A C ZP )\ 


dK m 


63 column: tkn [(£ e C$ Z + VC 2) Z ){C Z A C P 3 A C ZP R^ 

offers rV)+ 

M(c z k p c p ’ A c z ; A - c z ; a c p ; a c zp )\ 


(8.3.4-23d) 
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5u f'i row: 


pA column: 0 

<f>i r column: 0 

column: 0 

<j>i 3 column: 0 

UF"j column: 0 

9 F «j column: 0 


Uj column: 0 

u' g column: - tUm (C c C^ p ' + VCX' P ’ )C P ’ F ’ ^ 

h column: - 'k„UcC}? p ' + VC\'~ P ' )C PF ' 


(8.3.4-23e) 
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6r^F’ l, i row: 


jy 1 column: 0 

<j>f r column: 0 

<j)f 2 column: 0 

column: 0 


1 (r n W'Z' j_ T)r W'Z'\ r Z'F' r Z'F' 

u pn j column: f rrt ki{L c L lm + ^ 2 m ) u li 

(?/?-< y column: 0 

Uj column: c m kl(£cCu n Z + 'DC™ m z )C(< t C£j‘ 


tH ,( z\ r zr r z'F' 


u’ g column: c 3tt [( C + V< * )R% - M\C"C^ P ‘ 

+ vcZ z ')(c t Yc p Yc^ F U 

f-iZPf-iZF(~tP'F'\ nf 1 F' f^£ 

L/ pl L/ fco U ni f u F’o q u ii 


(3 


0 3 column: e mkp e ln3 (C c C^ z ' + VC\ v m z ' ){C ZP C Z J C p '/' 4 

c^ p c z fc p ; F ')R p F : r 


(8.3.4-23f) 


d*k 

du” a 
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Sui row: 


P/ 4 column: 0 

<})f r column: 0 

column: 0 

<j)y 3 column: 0 

UF n j column: 0 

9 ?" j column: 0 

Uj column: 0 

*i column: - + PC 1 T f " )C,? P ’ £= 

dup 

0 3 column: - t Mi {C c C^ P ' + PC 2 T P ' )C,f F ' 


(8.3.4-23g) 
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row: 


P / 1 column: 0 


row: 


column: 0 

<f)^ 2 column: 0 

0^3 column: 0 

upnj column: 0 

6p>fj column: 0 

ilj column: 0 

Up column: \{C c C\\ z 

0 3 column: 0 

column: 0 

<Pi r column: 0 

<^ 2 column: 0 

column: 0 

upnj column: 0 

Optij column: 0 

uj column: 0 

Up column: 0 




M )Cf/ 


d 2 «3 

du^du' 0 


(8.3.4-23h) 


(8.3.4-23i) 


0 3 column: 0 


The elements of E and F are determined from perturbations of Eqs. (8.3.4-15) which 
govern the lift, drag, and pitching moment. Thus the E matrix may be defined as: 


£ c row: 


Wz\ column: p a cciW Z i + \p* c ~t~^Z2 + 

/ CLQL 


XPaC 


, *2, ^ 
m 


T Z 12 


\V^ 2 column: p a cc x W^ 2 - + 


\W\ 212 


'KPaC 2 -7-^G < i 


G% 2 column: ^7rp 0 c 2 |W| 

£d 


(8.3.4-24a) 


T> row: 

Wz : i column: p a cc d W zl + \p aC ^^Z 2 

W% 2 column: p a cc d W$ 2 - \p aC ~^^zi 

Gzi 2 column: 0 

£ nc row: 

W z j column: 0 

Wz 2 column: 0 

G Q Z12 column: 0 


(8.3.4-24b) 


(8.3.4-24c) 
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M. row: 


W ® j column: p a c 2 c m Wz \ 4- -p«c 2 


i El ~q 
i« w iiVl Gzi2 


W§ 2 column: p a c 2 c m w£ 2 - - PaC 2 


. »'z. 


16 


7Tf) a C G 

ijyi 


Z12 


G^ 12 column: 
and for the elements of F 


1 3 

TTrt-r 

16 


W 3 |W| 


£ c row: 


P row: 


: Q 

W z j column: 0 


: <? 

11 ^ 2 column: 0 


G Z]2 column: 0 


Wz! column: 0 


W z 2 column: 0 


‘.Q 

G z \2 column: 0 


C nr row: 


: Q , irp a c 2 

W zi column: — - — 


: Q 

11 Z 2 column: 0 


G z 12 column: — 

zii 16 


(8.3.4-24d) 


(8.3.4-25a) 


(8.3.4-25b) 


(8.3.4-25c) 
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AA row: 


IV zi column: 


nPaC" 

16~ 


: Q 

W Z 2 column: 0 


(8.3.4-25d) 


>, g , 3 tt p a c 4 

G Z13 column: 


From the relations defining the relative velocity components and gradient, the elements of 
G and H can be determined. The elements of G are defined as 


low: 


Pf* column: 0 

cf)j r column: 0 

4>f 2 column: 0 

0^3 column: 0 


uf"> column: e a wC^ F ttF'ln C kj F ' ~ 


r 


+ cS*s^ A ) 


j3 


(8.3.4-26a) 


&F"j column: - fjktC F k F [C F A (U A 4- ry A ) + Vp!/j-t- 

n F F .',R%+ 


‘iu—t ??i A S?.f (cE A R% A + c£ a r% a )+ 

, _ riZ'F' n Z'F' c,F'I dP'F' 

f a tmfnojO w ll F , o n Fll 


Uj column: t ak iC^ m F fl F , T m C^ F ' 


^C F [ A (C% A R% A + C?, A R?, A ) 


F‘ A fjQ A\ 
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rtf r , 


u' column: - e Mm Cff + F^) + C^'V^)^ 


dK r 


C r P'F' r ZP 0 F'I nQP u *™ , 

Cal e 3/m^/fc ^33 Si F'fc^Z2 ^ 




+ c£^) 


<9*: r 


*»S 


n F’r&P , F' (r P'F' r ZP r Z’F’ , 
^ano^F'p^ JF' 7 V^/p ^ofc n</ ' 


yP'F' r>zPr< z 'F'\ ^ Km 

'lq C nfc ‘-'op 1 0 U " 


0 3 column: - CfcijC^f [C,^ '‘'(I/j' 4 + J’T' 4 ) + Cjn F ^F*n] + 


e w ,e a no(< F 'crcf/' + cg' , c*fcZ*')n£lRS:r 


fkl 3 


7 


/-tZPsiZ' A jjQP(/~iP' A nQ' A . /-> P' A nQ ^ 4 '\ 
C 2 i C a j ^Z2 V t/ *:2 -*^A2 + ^fc3 ir vl3 / 


G? 12 raw: 

Pj 4 column: 0 

< J)f r column: 0 

(^ 2 column: 0 

(/>^ 3 column: 0 


- t T y^f A f f ("'%* As~iF l A i /~iZ* A si F >l\ 

UF"j column: — C 13 [^ 22 u i 2 ^ ^23 °>3 / 

^{C^R% A + cfSn% J Hcf,'*s%* + c^RiS) 

r (8.3.4-26b) 
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'+ 


0F"j column: 

- A 

2 Z* A pQF r [/ (-iZ* AsyF'A i s~tZ* AsvF* A\ 

— l^ji n F , t |^o 22 o fc2 -j- ^ 23 is kz )- 

^ ffyZ* A pQ ‘A - fiZ' A pQ 1 A\/ fiF’ A pQ' A , pFM pQ ( >1\] 
Tyl C 22 K A2 + U 23 **,43 A°fc2 **yt2 + °fc2 U A3 /J 

- A 

7 f n Z'F' r F' A ir Z' A pQ' A , ^Z',4 pQ’A\ . 

"T" '° 22 i M2 + ° 23 ^3 / + 

r L 

+ c,r^)] } 


ij column: - ^Cf/' 1 [(Cf^C^ + Cg A c£ A )~ 

^(c£ A ii v M A + c^ A n^ A )(cf 2 ' A R% A + c£ a r* a )} 


-i 1 1-1 ZP^iP 1 F‘ n,F' l ^ Km 

u p column: €3 i m C 33 C lk U F , k 


m / A n PZ dQP \ (r ,Z‘A n P'A , r .Z'4 r *P'A\_ 

f tlm| t -'ll t - / /2 [1^22 '^Ic2 ' *-'23 *- / fc3 / 

+ c^ a r^ a )(c[; a r^ a + cz*#t*)] 
ch p cr; A {cZ A R<* 2 A + <£*</)+ 
c^ A c^{cr 2 ' A R^ A 4C,^<' 1 )]} 


r* 


s»? 


«j column: c,,,^ [(c£-‘c£ A + C?, A C& A )~ 

±(c£ A Ri 2 A + c£^xc£^ + 

/-» ZP f~*P* A ( (~iZ* A pQ* A . s-fZ* A nQ f v4\ 

°lfc ^11 V°22 **,42 + ^23 **,43 

cr, A c*f(cr; A Ri A + 


(8.3.4-26c) 
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The elements of H are defined as: 


W3 row: 


•„ A 

Pj column: 

n Z' A 

~ l Q 1 

lA 

0 Jr column: 

0 

lA 

0 12 column: 

r^Z' A f?Q A 

— C cri n A2 

lA 

</> 13 column: 

r Z'A joQ 1 A 
” C al K A3 

upnj column: 

r'Z'F' 

Opttj column: 

c nZ'F' nQP , , r Z'F' r Z'F' fiP'F' 

™Z2 ^fcm *^F r m 

Uj column: 

^rZ'F' 

u'p column: 

C r ZP nQP d* 3 
d “ J ° 33 du” 

column: 

C r ZPjjQP 
Z2 


(8.3.4-27a) 
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row: 


r>Q 

12 


■„ A 

P j column: 


<^ lr column: 
0 12 column: 
<^ 13 column: 


column: 


0/7^ column: 


“ G 11 G 22 


cf/c 


23 


0 

n Z'F' 
°3 j 


(8.3.4-27b) 


tZj column: 0 


u f 0 column: 


pZP 

U 33 


Qk-s 

du'> 


03 column: 

8.3.5. Spatial Discretization 

The variables U{ and #3 are expanded in a set of polynomials based on reference 30, 
The “CO” functions (u 3 and 63 ) are expanded in terms of x^i(x) where x = x$/L The 
functions used beyond the first two standard linear functions are orthonormalized. The 
Cl functions (u tt ) are expanded in terms of fli(x). The functions used beyond the first 
four standard cubic functions are orthonormalized. The details of the orthonormalization 
procedure are specified below. 

The expansions are given by 

t = l 


N 9 

u 3 =^2q 3i ^i(x) (8. 3. 5-1) 

t=l 

N A 

t=l 
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The functions V»i for i > 2 and 0i for i > 4 are constructed from the Jacobi polynomials 
g n (x) = G n _!(p,g;x) where p = 5 and g = 3 for the CO functions and where p = 9 and 
g = 5 for the Cl functions. 

Letting x = the CO shape functions are 

=1 — x 

ij) 2 —x (8. 3. 5-2) 

V’t — ®(1 - *)5«-2(®)/t-2 

where 3 < i < N + 1 and JV = JV 3 or iV 4 . The recursion relations used to compute the 
polynomials are 
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9\{x) = 1 


g 2 (x) =x - ^ 

9j{ x ) =92{x)gj-i{x) - gj- 2 (x)Aj-2 


5 ;(x) =0 


92 ( x ) =1 

9j( x ) =9j-i{ x ) + 92 {x)g' j ^ 1 {x) - g' j _ 2 (x)Aj- 2 


9 i( x ) =0 

92 ( x ) 

g'j{x) =2g' j _ 1 {x) + g 2 {x)g'^ 1 {x) - 9j_2 A j-2 

where 

t(i + 9 - l)(i+p-l){i+p-q) . * 

jfxj — <1 j r — V 

( 2 i 4 - p — 2 )( 2 i + p — l/( 2 i + p) 

/. 2 =30 



The derivatives of the shape functions are then 

v>; =-i 
+” =0 
t /4 =1 
V-2 =0 

V>I =[(1 - 2x)<7i_2 + *(1 - *)5i- 2 ]/i-2 

=(-25<_ 2 (x) - 4 (x - ^)flf'_ 2 (*) + *(! - x )g"-2( x )}fi-2 


(8. 3. 5-3) 


(8. 3. 5-4) 


(8. 3. 5-5) 
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Similarly, the Cl shape functions are 

/?j —1 — 3x 2 -\- 2x 

(3 2 =* - 2x 2 + x 3 

fa =3x 2 - 2x 3 (8. 3.5-6) 

fa =x 3 - i 2 


0, -x 2 (l - x) 2 g t - 4 (x)fi -4 

where 5 < i < N + 1 and JV = N } or N 2 . The g' s are the same as above for the CO shape 
functions, Ai is evaluated for p = 9, q = 5, and f 2 = 630. The higher derivatives are 

9i'( x ) =92 ) = 0 

g'!'{x) =3 + 92i x )9j-i ~ 9j'-2( x ) A i-2 


g""{x) =4 - 9j-2{ x ) A }-i 


and derivatives of the shape functions are 


0\ - — 6*(1 — x) 
fa = 3x 2 - 4x + 1 
0s = 6x(l - x) 
j3 f 4 — 3x 2 — 2;r 


/?" = 12(*-±) 

02 = 6 ( X - f) 

03 = -12(® - I ) 

04 = 6 ( X - I) 


/?;» = 12 = o 

p» = 6 = o 

/?'" = -12 = 0 

/?"' = 6 /C = 0 


(8. 3. 5-8) 


/?{ =[ — 4x(l — x)(x — + £ 2 (1 — ®) Sh- 4 ( x )]/«-4 

= [12(x 2 — ® + |)<7i_4(:r) - 8x(l — x)(x — |)<?I-4 + a;2 ( 1 “ x ) 3i-4( x )]/«-4 

/?'" = [24(* - \ )gi- a{x) + 36(x 2 - x + |)^- 4 ( x ) - 12x (l - x )( x - i)</i- 4 ( x )+ 

* 2 (1 - x) 2 g':U x ))fi-4 

0i" =j24^i_4(x) + 96(x — !)</i_4( a; ) + 72(x — x + g)<7j_ 4 ( a: ) — 

16x(l - x )( x - |)y''- 4 ( x ) + - x ) 2 9"- 4 ( x )]/i-4 

(8.3. 5-9) 

These formulas for shape functions, when substituted into expressions for virtual work 
of either internal, inertial, or applied loads, produce integrands that depend only on x 3 . 
These integrals can be evaluated to any accuracy desired by use of Gaussian quadrature. 
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8.3.6. Transformation from Nodal Coordinates 

In GRASP, a different set of generalized coordinates are used for the beam element 
than those for the nodes. It is therefore necessary to calculate the beam generalized 
coordinates in terms of the nodal displacements and rotational variables at both the root 
and tip of the beam, so that the beam equations can be written using a convenient set of 
generalized coordinates. 

The beam generalized coordinates q a i for i = 1,2, 3,4 determine the u Q displacements 
at the beam root and tip. Similarly, qt a determines the displacement at the root and 
tip, and q 4a determines the 6 3 rotation at the root and tip. The exact relations are 

qn ~um 


9l2 

r R'R 
— °31 

914 

“^31 

922 

y^rH'jR 
— ^32 

9 24 

r^T'R 

^32 

913 ^ 


923 > 

=C Uj 1 

932 i 



(8.3.6-1) 


</4 1 — sin 


-1 


q 4 2 — sin 


-1 





C T'R —(jT'TqTR 


The rotation expressions are derived from expressions for C p F written in terms of 
and #3 (see ref. 24, equations 4, 17, and 60-62) for which 

¥ 3 # I ¥71 1 # . 

C/j r = sin &2 = Uj 


p # f rpH t | 

C 32 — cos 62 sm#i = u 2 

Cfj F = — cos &2 sin 03 — — \J 1 — u'j 2 sin 63 


(8. 3. 6-2) 
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8.3.7. Transformation to Forces and Moments 

The generalized forces calculated for the beam element root and tip correspond to 
the beam generalized coordinates. These forces must now be transformed into forces and 
moments at the root and tip nodes. The virtual work at the root can be written in terms of 
the static residuals Qr and the linear coefficient matrices. In terms of the beam generalized 
coordinates, this relation is 


-<5Wfl = 6qR T { — Qr + L R qn) 


(8.3. 7-1) 


where 8q R T = [^9n ^921 ^931 ^912 ^922 ^94iJ an d 9 r — l_9ii 921 93i 9i2 922 94iJ> aI1( l 
is a linear operator representing {Mr ^r)- Note that this equation defines 

the negative of the virtual work. The explanation for treating the virtual work in this 
manner is that it is conventional for Lr to be positive, and Lr is normally considered to 
be positive on the left-hand side of the equations of motion, while Qr is positive on the 
right-hand side. 

The root node virtual displacements and rotations may be related to the beam virtual 
generalized coordinates by the expression 


8qu ' 


r 



£921 


A 

0 

f Sur 

8qs 1 

> = 

3x3 

3x3 

3X! 

8q 12 

0 

Rr 

] Hr 

8q 22 


3x3 

3x3 

l 3X1 J 

6941 j 






(8.3. 7-2) 


where the root node virtual displacements are 8ur = R , and the root node virtual 

rotations are 8ip R — Hr R - The 6x6 coefficient matrix that premultiplies the root node 
virtual displacement and rotation vector is called Tr and matrix Rr (ref. 24, eq. 67) is 
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Similarly, the perturbed root node displacements ur and rotations Or are related to 
the perturbed element generalized coordinates q R through the expression 


Or 



(8. 3. 7-4) 
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When the virtual work at the beam root is transformed into nodal coordinates by 
the substitution of equations (8. 3. 7-2) and (8. 3. 7-4) into equation (8. 3. 7-1), the following 
expression is obtained: 


-6Wr — [6ur^ 
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(8. 3. 7-5) 


First consider the transformation of Lr, which contains the dynamic matrices Mr, 
Cr , and Kr. The transformation of the element generalized coordinates into the nodal 
generalized coordinates introduces the transformation matrix Tr into the expression for the 
virtual work. Since Tr is a function of C R which is a function of the nodal rotations, 
it must also be perturbed to recover any additional perturbation contributions. In the 
case of the linearized dynamic matrices A//?, Cr , and Kr, no new perturbation terms are 
introduced by the transformation, since any such contributions would be nonlinear. 

The transformation of the static generalized force Qr is, however, another matter. In 
this case, transformation does contribute an additional term, called the geometric stiffness 
term K to the linearized dynamic equations. Geometric stiffness originates from the 
perturbation of Tr. 

-6q rQr — -[Sur T S^r T \Tr T {Qr} (8. 3. 7-6) 

When Tr is perturbed, 
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When equation (8. 3. 7-7) is multiplied out, only one of the 3x3 submatrices is nonzero. 
This submatrix is called the root geometric stiffness matrix kcR } and it contains only terms 
that are related to the nodal rotations. 
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where 


(8.3.7-10) 


Qr = 



Written in index notation to allow the isolation of k G R, the root geometric stiffness is 

T~ -.ijRiQne* (8.3.7-11) 
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The geometric stiffness matrix used to transform all of the root nodal degrees of freedom 
is then 
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The virtual work at the root can now be written in the form 
-6Wr = [Sur T 6i/>r T J 1-Tr T Qr + T r 1 LrTr |^| + Kr (-J} 

= [Sur T bl^R 1 J Qr + L*r | | 

where _ 

Qr =Tr T Qr 

L*r =Tr t LrTr + Tr t k% 


(8.3.7-14) 


(8.3.7-15) 


The transformation of the generalized forces and moments at the tip of the element 
into nodal forces and moments is similar to that for the root. In beam element generalized 
coordinates, the virtual work at the tip is 

— SWt = $<1t T {~Qt d- (8.3.7-16) 


where SqT^ — |^9i3 ^923 ^932 ^9i4 ^924 ^942] an d <7^ — [913 923 932 9i4 924 942 Ji 
Lt is a linear operator representing — + K T ). Note the similarity with 

equation (8. 3. 7-1). 
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The equation that relates the tip node virtual displacements to the element virtual 
generalized coordinates is 


(8.3.7-17) 


where the tip node virtual displacements are bur — bR^ n , and the tip node virtual 
rotations ore btfrr = b^]^ R . The 6x6 coefficient matrix that premultiplies the tip node 
virtual displacement and rotation vector is called TV and matrix Rt is 
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Similarly, the perturbed tip node displacements ur and rotations Or are related to 
the perturbed element generalized coordinates qr through the expression 
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The expression for the virtual work at the tip is similar to the expression for the 
virtual work at the root. 

-SW T = 1 &u t T Ht T J {{-T t t Q t } + T t t L t Tt {]£} } (8.3.7-20) 


As in the derivation of the transformation of Lr , no additional terms result from the 
transformation of Lr- There are, however, geometric stiffness terms that result from the 
transformation of Qr- Following the derivation of the root geometric stiffness, 
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When equation (8.3.7-21) is expanded, only one, nonzero 3x3 submatrix remains. It is 
called the tip geometric stiffness matrix k G r , and 
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Written 

in index notation to allow the isolation of fccT> the tip geometric 

stiffness is 
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The virtual work at the tip can now be written in the form 
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For both the root and tip, derivatives of Rr and R T with respect to the q i} are needed. 
The only nonzero elements of these arrays may determined from 
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C nR at the root and C T n at the tip and R is Rr at the root and Rx at the 
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9. CONCLUDING REMARKS 


In response to the limitations of previous methods for analyzing rotorcraft, GRASP 
has been developed. GR ASP is a general-purpose program which treats the nonlinear static 
and linearized dynamic behavior of rotorcraft represented by arbitrarily connected rigid- 
body and beam elements. Large relative motions and deformation-induced displacements 
and rotations are permitted (as long as the strains in the beam element are small). Periodic 
coefficients are not treated, restricting the solutions to rotorcraft in axial flight and on the 
ground. 

GRASP uses a modern approach for modeling structures, incorporating the features 
of several traditional methods. The basic approach which provides the foundation for 
large relative motion kinematics is derived from “multibody” research with an expanded 
emphasis on multiple levels of substructures. This is combined with the finite element 
approach which provides flexible modeling through the use of libraries of elements, con- 
straints, and nodes. The use of a variable-order polynomial beam element makes the finite 
element approach more effective. The incorporation of aeroelastic effects, including inflow 
dynamics and nonlinear aerodynamic coefficients for the beam element, further extends 
the capabilities. 

Due to the fact that GRASP was developed using structured, modular, software meth- 
ods, changes to the code are relatively easy to perform. This makes it practical to modify 
the code in order to enhance its functionality. Some of the many areas where possibilities 
for enhancements exist are expanded solution procedures, improved aerodynamic models, 
expanded modeling capabilities, new elements, and new constraints. 

Existing solution procedures (steady-state and asymmetric eigenproblem) could easily 
be expanded to include a symmetric eigensolution. This solution procedure would take the 
symmetric part of the linearized, perturbed equations of motion, then calculate the eigen- 
values and eigenvectors. The symmetric eigensolution would be to generate the modes for 
another new solution procedure, the subspace reduction. The subspace reduction would 
allow the user to solve for the asymmetric eigensolution using a reduced set of admissible 
functions. A reference deformations solution procedure would allow a user to take any 
steady-state solution and use it either as an initial guess for another steady-state problem, 
or as the state about which the linearization is performed for an eigensolution. The refer- 
ence deformations solution would lift the restriction that the same model must be used in 
the the steady-state solution and the eigensolution. Another valuable enhancement would 
be to extend GRASP to forward flight using either a time-domain solution, a periodic 
solution, or both. 

Enhancements to the aerodynamics could include adding the capability for table- 
lookup for the aerodynamic coefficients, and perhaps making those coefficients functions 
of Mach number. Another possibility would be to incorporate a lifting-line or lifting surface 
theory to calculate the aerodynamic forces. Wake geometry could also be included. Other 
valuable enhancements to the aerodynamic, model would be the inclusion of transonic and 
dynamic stall effects. 
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The modeling capabilities could also be improved with the addition of the ability to 
model applied loads. It might be advantageous to include simple, dead loads (forces), 
and geometrically nonlinear loads such as applied moments and nonconservative forces. 
With the rapid growth of control theory, some sort of control representation should be 
included in GRASP. This could be as simple as specifying the thrust level of the rotor, or 
as complex as a complete control representation including sensors, actuators, and control 
laws. In addition, it would be convenient to implement a “generic” node. Such a node 
would be used to allow the user to define generalized coordinates not associated with any 
of the predefined nodes. 

GRASP would greatly benefit from the addition of a composite beam element and 
a direct-input element. The composite beam element would be able to rigorously treat 
the structural couplings introduced by composite layups. This element might also include 
the effects of shear deformation, initial curvature, and warping rigidity. The direct-input 
element would be used in conjunction with the generic node to allow the user to define the 
properties of elements that are not included in GRASP. An example of such a use would be 
taking a set of modes from a NASTRAN analysis to represent the fuselage of a helicopter. 

New constraints that would enhance the capabilities of GRASP include a moving- 
frame constraint, a pin constraint, and a clamp constraint. The moving-frame constraint 
would allow a frame to deform with the structure. Currently, frame motion is independent 
of the structure. The pin constraint would allow a node to rotate arbitrarily about either 
a frame or another node. Eliminating all motion of a node would be accomplished using 
the clamp constraint. 

From this description of possible enhancements, it should be obvious that GRASP 
has a great potential for growth. Because of its modular construction, GRASP has the 
capability to handle expansion without requiring massive rewriting of the existing equations 
and code. This framework makes GRASP a desirable platform for future development. 
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